{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os\n",
    "import re\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "from scipy.stats import spearmanr, pearsonr\n",
    "from scipy.stats import t, norm\n",
    "from tqdm import tqdm\n",
    "\n",
    "\n",
    "def get_time(x):\n",
    "    if x == 'UT':\n",
    "        return x\n",
    "    else:\n",
    "        pattern = re.compile(r'\\d+h')\n",
    "        return re.findall(pattern, x)[0]\n",
    "\n",
    "\n",
    "class DATASET:\n",
    "    def __init__(self, datasetname):\n",
    "        self.name = datasetname\n",
    "        self.path_prefix = Path(\"./seurat_objects\")\n",
    "        self.information = self.get_information()\n",
    "    def get_information(self):\n",
    "        if self.name == 'onemillionv2':\n",
    "            self.path = '1M_v2_mediumQC_ctd_rnanormed_demuxids_20201029.sct.h5ad'\n",
    "            self.individual_id_col = 'assignment'\n",
    "            self.timepoint_id_col = 'time'\n",
    "            self.celltype_id = 'cell_type_lowerres'\n",
    "            self.chosen_condition = {'UT': 'UT',\n",
    "                                     'stimulated': '3h'}\n",
    "        elif self.name == 'onemillionv3':\n",
    "            self.path = '1M_v3_mediumQC_ctd_rnanormed_demuxids_20201106.SCT.h5ad'\n",
    "            self.individual_id_col = 'assignment'\n",
    "            self.timepoint_id_col = 'time'\n",
    "            self.celltype_id = 'cell_type_lowerres'\n",
    "            self.chosen_condition = {'UT': 'UT',\n",
    "                                     'stimulated': '3h'}\n",
    "        elif self.name == 'stemiv2':\n",
    "            self.path = 'cardio.integrated.20210301.stemiv2.h5ad'\n",
    "            self.individual_id_col = 'assignment.final'\n",
    "            self.timepoint_id_col = 'timepoint.final'\n",
    "            self.celltype_id = 'cell_type_lowerres'\n",
    "            self.chosen_condition = {'UT': 't8w',\n",
    "                                     'stimulated': 'Baseline'}\n",
    "        elif self.name == 'ng':\n",
    "            self.path = 'pilot3_seurat3_200420_sct_azimuth.h5ad'\n",
    "            self.individual_id_col = 'snumber'\n",
    "            self.celltype_id = 'cell_type_mapped_to_onemillion'\n",
    "        else:\n",
    "            raise IOError(\"Dataset name not understood.\")\n",
    "    def load_dataset(self):\n",
    "        self.get_information()\n",
    "        print(f'Loading dataset {self.name} from {self.path_prefix} {self.path}')\n",
    "        self.data_sc = sc.read_h5ad(self.path_prefix / self.path)\n",
    "        if self.name.startswith('onemillion'):\n",
    "            self.data_sc.obs['time'] = [get_time(item) for item in self.data_sc.obs['timepoint']]\n",
    "        elif self.name == 'ng':\n",
    "            celltype_maping = {'CD4 T': 'CD4T', 'CD8 T': 'CD8T', 'Mono': 'monocyte', 'DC': 'DC', 'NK': 'NK',\n",
    "                               'other T': 'otherT', 'other': 'other', 'B': 'B'}\n",
    "            self.data_sc.obs['cell_type_mapped_to_onemillion'] = [celltype_maping.get(name) for name in\n",
    "                                                          self.data_sc.obs['predicted.celltype.l1']]\n",
    "    def get_cMono_ncMono(self):\n",
    "        def tell_cmono_foronemillion(x):\n",
    "            if x == 'mono 1' or x == 'mono 3' or x == 'mono 4':\n",
    "                return 'cMono'\n",
    "            elif x == 'mono 2':\n",
    "                return 'ncMono'\n",
    "        if self.name.startswith('onemillion'):\n",
    "            self.data_sc.obs['sub_monocytes'] = [tell_cmono_foronemillion(x) for x in\n",
    "                                                 self.data_sc.obs['cell_type']]\n",
    "            self.cmono = self.data_sc[self.data_sc.obs['sub_monocytes'] == 'cMono']\n",
    "            self.ncmono = self.data_sc[self.data_sc.obs['sub_monocytes'] == 'ncMono']\n",
    "        elif self.name.startswith('stemi'):\n",
    "            self.cmono = self.data_sc[self.data_sc.obs['cell_type'] == 'cMono']\n",
    "            self.ncmono = self.data_sc[self.data_sc.obs['cell_type'] == 'ncMono']\n",
    "        elif self.name == 'ng':\n",
    "            self.cmono = self.data_sc[self.data_sc.obs['predicted.celltype.l2'] == 'CD14 Mono']\n",
    "            self.ncmono = self.data_sc[self.data_sc.obs['predicted.celltype.l2'] == 'CD16 Mono']\n",
    "        else:\n",
    "            raise IOError(\"Dataset name not understood.\")\n",
    "\n",
    "example_savedir = Path(\n",
    "    \"/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/output/examples\"\n",
    ")\n",
    "\n",
    "import subprocess\n",
    "bashfile_path = '/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/bios/select_snps_from_vcf.sh'\n",
    "def get_snps_from_vcffile(bashfile_path, vcf_path, snps_path, savepath):\n",
    "    response = subprocess.run([bashfile_path, vcf_path, snps_path, savepath])\n",
    "    print(response)\n",
    "    return None\n",
    "\n",
    "# sample id mapping\n",
    "gtefile = pd.read_csv(\n",
    "    '/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/input/summary/gte-fix.tsv',\n",
    "    sep='\\t'\n",
    ")\n",
    "gte_dict = gtefile.set_index(\"expressionsampleID\")[\"genotypesampleID\"].T.to_dict()\n",
    "\n",
    "\n",
    "def corr_to_z(coef, num):\n",
    "    t_statistic = coef * np.sqrt((num - 2) / (1 - coef ** 2))\n",
    "    prob = t.cdf(t_statistic, num - 2)\n",
    "    z_score = norm.ppf(prob)\n",
    "    positive_coef_probs = 1 - prob\n",
    "    positive_coef_probs[coef < 0] = 0\n",
    "    negative_coef_probs = prob\n",
    "    negative_coef_probs[coef > 0] = 0\n",
    "    probs = negative_coef_probs + positive_coef_probs\n",
    "    return z_score, probs\n",
    "\n",
    "\n",
    "def get_individual_networks_selected_genepairs(data_df, data_sc, individual_colname, genepair, fillna=False):\n",
    "#     data_df = pd.DataFrame(data=data_sc.X.toarray(),\n",
    "#                            index=data_sc.obs.index,\n",
    "#                            columns=data_sc.var.index)\n",
    "    gene1, gene2 = genepair.split(';')\n",
    "    sorted_genepair = [';'.join(sorted([gene1, gene2]))]\n",
    "    coef_df = pd.DataFrame(index=sorted_genepair)\n",
    "    coef_p_df = pd.DataFrame(index=sorted_genepair)\n",
    "    zscore_df = pd.DataFrame(index=sorted_genepair)\n",
    "    zscore_p_df = pd.DataFrame(index=sorted_genepair)\n",
    "    data_selected_df = data_df[[gene1, gene2]]\n",
    "    print(\n",
    "        f\"Calculating networks for {len(data_sc.obs[individual_colname].unique())} individuals and;\\n{genepair}\"\n",
    "    )\n",
    "    for ind_id in tqdm(data_sc.obs[individual_colname].unique()):\n",
    "        cell_num = data_sc.obs[data_sc.obs[individual_colname] == ind_id].shape[0]\n",
    "        if cell_num > 10:\n",
    "            individual_df = data_selected_df.loc[data_sc.obs[individual_colname] == ind_id]\n",
    "            individual_coefs, individual_coef_ps = spearmanr(individual_df.values, axis=0)\n",
    "            if data_selected_df.shape[1] == 2:\n",
    "                individual_coefs_flatten = pd.DataFrame(data = [individual_coefs],\n",
    "                                                        index = sorted_genepair)\n",
    "                individual_coef_ps_flatten = \\\n",
    "                pd.DataFrame(data=[individual_coef_ps],\n",
    "                             index=sorted_genepair)\n",
    "            else:\n",
    "                individual_coefs_flatten = pd.DataFrame(\n",
    "                    data=individual_coefs[np.triu_indices_from(individual_coefs, 1)],\n",
    "                    index=sorted_genepair).loc[sorted_genepair]\n",
    "                individual_coef_ps_flatten = \\\n",
    "                pd.DataFrame(data=individual_coef_ps[np.triu_indices_from(individual_coefs, 1)],\n",
    "                             index=sorted_genepair).loc[sorted_genepair]\n",
    "            coef_df[ind_id] = individual_coefs_flatten\n",
    "            coef_p_df[ind_id] = individual_coef_ps_flatten\n",
    "            try:\n",
    "                individual_zscores_flatten, individual_zscore_ps_flatten = corr_to_z(\n",
    "                    individual_coefs_flatten.values,\n",
    "                    cell_num\n",
    "                )\n",
    "                zscore_df[ind_id] = individual_zscores_flatten\n",
    "                zscore_p_df[ind_id] = individual_zscore_ps_flatten\n",
    "            except:\n",
    "                continue\n",
    "        else:\n",
    "            print(\"Deleted this individual because of low cell number\", cell_num)\n",
    "    if fillna:\n",
    "        zscore_df = zscore_df.fillna(0)\n",
    "    return data_selected_df, zscore_df, zscore_p_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading dataset onemillionv2 from /groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/seurat_objects 1M_v2_mediumQC_ctd_rnanormed_demuxids_20201029.sct.h5ad\n"
     ]
    }
   ],
   "source": [
    "datasetname = 'onemillionv2'\n",
    "dataset = DATASET(datasetname)\n",
    "dataset.load_dataset()\n",
    "data_sc = dataset.data_sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CompletedProcess(args=['/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/bios/select_snps_from_vcf.sh', '/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/output/genotypevcfs/chr1/GenotypeData.vcf.gz', PosixPath('/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/output/examples/snplist.rs221045'), PosixPath('/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/output/examples/rs221045.vcf')], returncode=0)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>#CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>ID</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th>QUAL</th>\n",
       "      <th>FILTER</th>\n",
       "      <th>INFO</th>\n",
       "      <th>FORMAT</th>\n",
       "      <th>1_LLDeep_1191</th>\n",
       "      <th>...</th>\n",
       "      <th>s21</th>\n",
       "      <th>s43</th>\n",
       "      <th>s24</th>\n",
       "      <th>s23</th>\n",
       "      <th>s45</th>\n",
       "      <th>s26</th>\n",
       "      <th>s25</th>\n",
       "      <th>s28</th>\n",
       "      <th>s27</th>\n",
       "      <th>s29</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>16530049</td>\n",
       "      <td>rs221045</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>GT:DS</td>\n",
       "      <td>0/0:0.03</td>\n",
       "      <td>...</td>\n",
       "      <td>0/1:1.0</td>\n",
       "      <td>0/0:0.010000000000000009</td>\n",
       "      <td>0/1:1.0</td>\n",
       "      <td>0/0:0.0</td>\n",
       "      <td>0/0:0.0</td>\n",
       "      <td>1/1:2.0</td>\n",
       "      <td>0/0:0.0</td>\n",
       "      <td>0/1:1.0</td>\n",
       "      <td>0/0:0.0</td>\n",
       "      <td>0/1:1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 182 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   #CHROM       POS        ID REF ALT QUAL FILTER INFO FORMAT 1_LLDeep_1191  \\\n",
       "0       1  16530049  rs221045   T   C    .      .    .  GT:DS      0/0:0.03   \n",
       "\n",
       "   ...      s21                       s43      s24      s23      s45      s26  \\\n",
       "0  ...  0/1:1.0  0/0:0.010000000000000009  0/1:1.0  0/0:0.0  0/0:0.0  1/1:2.0   \n",
       "\n",
       "       s25      s28      s27      s29  \n",
       "0  0/0:0.0  0/1:1.0  0/0:0.0  0/1:1.0  \n",
       "\n",
       "[1 rows x 182 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "celltype = 'monocyte'\n",
    "snp_id = 'rs221045'\n",
    "chromosome = '1'\n",
    "snp_vcf_path = example_savedir/f'{snp_id}.vcf'\n",
    "with open(example_savedir/f'snplist.{snp_id}', 'w') as f:\n",
    "    f.write(f'{snp_id}\\n')\n",
    "vcf_path = f'/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/output/genotypevcfs/chr{chromosome}/GenotypeData.vcf.gz'\n",
    "get_snps_from_vcffile(bashfile_path, vcf_path, example_savedir/f'snplist.{snp_id}', snp_vcf_path)\n",
    "gt = pd.read_csv(snp_vcf_path, sep='\\t', skiprows=6)\n",
    "gt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating networks for 72 individuals and;\n",
      "AC005076.5;ARHGEF19\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/72 [00:00<?, ?it/s]/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/tools/Beeline/miniconda/envs/scpy3.8/lib/python3.8/site-packages/scipy/stats/stats.py:4264: SpearmanRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.\n",
      "  warnings.warn(SpearmanRConstantInputWarning())\n",
      "100%|██████████| 72/72 [00:00<00:00, 210.51it/s]\n",
      "/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/tools/Beeline/miniconda/envs/scpy3.8/lib/python3.8/site-packages/seaborn/categorical.py:1296: UserWarning: 42.5% of the points cannot be placed; you may want to decrease the size of the markers or use stripplot.\n",
      "  warnings.warn(msg, UserWarning)\n",
      "/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/tools/Beeline/miniconda/envs/scpy3.8/lib/python3.8/site-packages/seaborn/categorical.py:1296: UserWarning: 7.1% of the points cannot be placed; you may want to decrease the size of the markers or use stripplot.\n",
      "  warnings.warn(msg, UserWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Not Imputed SpearmanrResult(correlation=-0.028018282506059713, pvalue=0.8942369051146191)\n",
      "Imputed SpearmanrResult(correlation=-0.24638574744096847, pvalue=0.03833253459364005)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAFNCAYAAABbpPhvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACfvklEQVR4nOydd3hTZdvAfyez6d67lNFS9kY2KAgoQ8SJ+vm+6Kso4sa99x6v8oqICu49QERkCgjKBtl7QxfdTdvM8/2RNm1o0zZtkybt87suLpIznufOac6d+9zPPSRZlmUEAoFAIBAIBF6ForkFEAgEAoFAIBBURxhpAoFAIBAIBF6IMNIEAoFAIBAIvBBhpAkEAoFAIBB4IcJIEwgEAoFAIPBChJEmEAgEAoFA4IUII00AQFpaGidOnGhuMVoEjzzyCG+//XZziyEQCJywceNGhg8f3txitBhGjhzJX3/91dxitEiEkdaENNcX9aeffuK6667z+LwNZdGiRVx00UX06tWLO+64g/z8/BqPy8nJ4f7772fo0KH07duXKVOm8M8//9j3b9y4kU6dOtG7d2/7v59//tlDn6J5qO+1Azh9+jQ33ngjPXv25JJLLnH4bm7YsIGJEyfSr18/BgwYwIwZM8jMzLTvz8/P595772XAgAEMGDCAmTNnUlxcDEBubi5TpkxhwIAB9OvXj2uvvZatW7e67TMLGsbIkSMZPHgwJSUl9m3ff/89N954Y73Ov/HGG/n++++d7j99+jRpaWmYzeZGy+oqvvQgJMsyr7/+uv1eeu2116itPOnff//NJZdcQs+ePbnxxhs5c+aMfd8nn3zCqFGj6NOnD0OHDuWll15qluvvKVy9dr/99huXXnopvXv3Zty4caxYscK+r7CwkIcffphBgwYxaNAgZs2aVeMYmzZtIi0tzeH7tXr1aq677jr69evHkCFDeOKJJ+z60N0II03gErIsY7VaG3z+oUOHeOqpp3jttddYv349Op2OZ599tsZjS0pK6N69Oz/99BObNm1i8uTJTJs2Db1ebz8mOjqa7du32/9Nnjy5wbJ5O65cO4CZM2fSpUsXNm7cyH333cfdd99Nbm4uACkpKXz00Uds2bKFP//8k+TkZJ5++mn7uf/9738pLCxk5cqVrFixgpycHLtSCwgI4KWXXuLvv/9m8+bN3HrrrUyfPr1F/1j4KhaLhc8++6y5xfBpGvu9/vbbb1mxYgULFy7kl19+YfXq1XzzzTc1Hpubm8udd97JPffcw6ZNm+jWrRv33Xefff/IkSP5+eef2bZtG7/++iv79+/n888/b5R83owr1y4zM5OHHnqIRx55hG3btvHQQw8xc+ZMcnJyAHj55ZcpLS1l1apVfP/99yxcuJAff/zRYQyTycSLL75Iz549HbYXFRUxffp0/vzzT3777TcyMjJ47bXX3POhz0MYaW7ip59+YsqUKbz00kv069ePUaNGsW3bNn766SdGjBjBoEGDHLw+jzzyCE899RQ33XQTvXv35v/+7//sT1A1PbFWPOUeOXKEp59+mh07dtC7d2/69esHgNFo5NVXX+XCCy9k8ODBPPXUU5SVldnP/+ijjxg6dChDhw7lhx9+qPWz3Hjjjbz99ttMmTKFnj17curUqQZfl0WLFjFy5Ej69+9PQEAA99xzD8uXL6/xqSQpKYmbbrqJ6OholEol1157LSaTiWPHjtVrrrlz53Lbbbc53T9y5Eg++OADxo0bR//+/Xn00UcxGAwAXHrppfzxxx/2Y81mMwMGDGDPnj0A3H333QwZMoS+fftyww03cOjQoRrnqMnLWXVpua6/U1VcuXbHjh1jz5493HXXXfj5+TF27Fg6duzI0qVLAYiMjCQmJsZ+vFKp5OTJk/b3p0+fZtSoUQQGBhIUFMTo0aM5fPgwAFqtlvbt26NQKJBlGYVCQUFBAQUFBU6vtaB5+M9//sO8efMoLCyscf+2bdu48sor6du3L1deeSXbtm0D4O2332bLli0899xz9O7dm+eee67OuR555BGeeeYZbrnlFnr37s2UKVPIzs7mxRdfpH///lxyySXs3bvXfnxt919t9823337LokWL+Pjjj+nduze33347YPuRvuuuuxg4cCAjR450ME7Lysp45JFH6N+/P+PGjWPXrl21fpa0tDS+/PJLxowZw5gxY+r87LWxYMECbr75ZmJjY4mJieGmm25y6vFfvnw5qampXHrppWi1Wu666y7279/PkSNHAGjTpg3BwcEA9nuvapjKbbfdxty5c2scu+J35Ntvv7Xr/nnz5gG2a9ejRw8Hz/zevXsZMGAAJpOJkydP8q9//cvBs+7sO3W+l/P8peXa/k6NuXYZGRkEBQUxYsQIJEniwgsvRKfT2fXaqlWruOWWW9DpdCQmJnLVVVdVM9Lmz5/PkCFDaN++vcP2iRMnMnz4cHQ6HSEhIVxzzTVs377dqdxNiTDS3MjOnTtJS0tj48aNTJgwgfvvv59du3axfPlyXn/9dZ577jkHr9CiRYu444477Mt4DzzwQJ1zdOjQgWeffZZevXqxfft2tmzZAsDrr7/OsWPHWLBgAcuWLSMrK4v33nsPgLVr1zJv3jzmzZvHsmXL+Pvvv+ucZ+HChTz//PNs27aN+Ph4nnnmGfr161fjv4kTJzod59ChQ6Slpdnft2nTBrVazfHjx+uUYd++fZhMJpKTk+3bcnNzGTx4MCNHjuSll15yWNqZNm0aH3zwQa1jVij75cuXc+zYMWbPng3A+PHj+fXXX+3HrVu3jrCwMLp27QrA8OHDWbp0KX///TddunSp19+qJmr7O52PK9fu8OHDJCUlERgYaN/WqVMnu6EFcPbsWfr160ePHj2YN28et9xyi33fDTfcwOrVq+3G19KlSxk2bJjDHBMnTqRHjx5Mnz6dq6++moiIiAZdA4H76NatGxdccAEff/xxtX35+fncdttt3HjjjWzcuJGbbrqJ2267jby8PO677z769evHU089xfbt23nqqafqNd+SJUu499572bBhAxqNhmuvvZauXbuyYcMGxo4dy8svv+xwvLP7rzauvfZaJk6cyH/+8x+2b9/OnDlzsFqtTJ8+nbS0NNauXcunn37Kp59+yp9//gnA//73P06ePMny5cv5+OOPWbBgQZ3zrFixgu+++47ffvsNwB4eUNO/Z555xuk4hw4dolOnTvb3nTp1cvpQd/497u/vT5s2bRzu20WLFtGnTx8GDhzI/v37mTJlin3fBx98wLRp02r9XBs3bmTZsmV8/PHHzJ07l7/++ouYmBh69erFsmXLHOYZO3YsarUaWZa57bbb+PPPP1myZAkZGRlOlwtro66/U03Xo77Xrlu3bnTo0IGVK1disVhYsWIFGo3G4XpWRZZlh7HOnDnDjz/+yIwZM+r8HJs3byYlJaXO45oCYaS5kcTERK688kqUSiXjxo0jPT2dGTNmoNFoGDp0KBqNxsF7ceGFF9K/f380Gg333XcfO3bsID093eV5ZVnm+++/57HHHiM0NJTAwEBuu+02Fi9eDNgU6RVXXEHHjh3x9/fnzjvvrHPMyZMnk5qaikqlQq1W88wzz7Bly5Ya/y1atMjpOCUlJQQFBTlsCwwMdDBWa6K4uJiHHnqIO++8035++/btWbBgAevWrePTTz9lz549vPLKK3V+lqrccMMNxMXFERoayvTp0+3XaOLEiaxatYrS0lLAprAmTJhgP++qq64iMDAQjUZjf9otKipyae66/k7n48q10+v11Y4NCgpyODY+Pp4tW7awYcMG7rnnHoenxy5dumAymexPzkqlkuuvv95hvEWLFrF161befPNN+vbt69JnF3iOu+++my+++MK+1F3B6tWrSU5O5vLLL0elUjFhwgTat2/v4EF2ldGjR9OtWze0Wi2jR49Gq9Vy+eWX23Xgvn37HI53dv+5yq5du+xLhRqNhqSkJK655hq7gbVkyRJuv/12QkNDiYuLq1dc3rRp0wgNDcXPzw+wfd+d6bzajLSSkhKHh6WgoCBKSkpqjK2qzz0+ceJEtm3bxtKlS5kyZYrLD0czZszA39+ftLQ0rrjiCvvD6MSJE+2vZVnmt99+sz9wJycnM2TIEDQaDeHh4dx0001s3rzZpXmh7r/T+bhy7ZRKJZMmTeKBBx6ge/fuzJw5k+eeew5/f38Ahg0bxty5cykuLubEiRP8+OOPdv0O8MILL3DPPfcQEBBQ62dYv349CxYs4O6773b58zcElUdmaaVUvXkqbvTIyEj7Nq1W63DzxcbG2l8HBAQQEhJCVlaWyzdhbm4upaWlXHHFFfZtVWPJsrKy6Natm31fQkJCnWPGxcW5JAPAli1buPXWWwGbQbB48WL8/f2rLc8VFxfXemOUlZVx++2307NnT4fly6ioKKKiogDb0uiDDz7IbbfdVq+lmQqqfq74+HiysrIAm1Lq0KEDf/zxBxdddBGrVq2yP31bLBbefvttfv/9d3Jzc1EobM86eXl51RRsbdT1dzofV65dQEBAvY8NDQ1l8uTJTJo0ibVr16JSqbjnnnvo1KkTs2fPRpZlXn31VR588EHeeecdh3O1Wi0TJkzg0ksvpXPnzg5PvQLvoGPHjlx44YXMnTuXDh062LdnZWURHx/vcGx8fLxDAomrnK/zquo7Pz8/B083OL//XOXMmTNkZWXZwz3Adp9WvM/Kyqo2V100ROfNmTPH7r2fOHGi3UioqueLi4vx9/dHkqRq59d0j+v1+hrv27Zt25Kamsqzzz7L//73v3rLWPVzJSQkcPDgQQDGjh3L888/T2ZmJidOnECSJPv1y8nJ4YUXXmDLli3o9XpkWbYvu7pCXX+n83Hl2v3111+88cYbfPbZZ3Tt2pXdu3dzxx138OGHH9K5c2eeeOIJnn/+ecaOHUtoaCjjx4+3PxSsWrUKvV7PuHHjapV/x44dzJw5k3fffZd27dq5/PkbgjDSvIiMjAz7a71eT0FBAdHR0Wi1WsBmrFQ8VWRnZ9uPPf8LGxYWhp+fH4sXL3aIO6ogOjrawUN39uzZOmU7f46nnnrKqceswiDr169ftXX71NRU9u/fb39/6tQpTCYTbdu2rXEso9HIjBkziImJqdP4kiSp1syfmjj/OkRHR9vfT5gwgV9//RWr1UpKSop9mXXRokWsXLmS+fPnk5iYSFFREf37969xbp1O5xBjVvXvVtff6XxcuXYpKSmcOnWK4uJi+3dm//79Dt7AqlgsFnJyciguLiY0NJQDBw7wzDPP2J9Cr7vuumqetKqYzWZOnToljDQv5e6772by5MncfPPN9m3R0dHV7v309PRqy9ruxNn9V9t9A9X1UVxcHImJiQ7LdVWJiooiPT2d1NTUavM64/w5xo8f71RXVhhkt99+uz1GroKK+7ZHjx6A7T6skON8UlNTHWKuSkpKOHnypNOlNbPZ7LAaUx/S09PtxnrVax4cHMyQIUNYsmQJR48eZfz48fZr8OabbyJJEr/88gthYWGsWLHCqT4+/2937tw5++u6/k7n48q127dvH/369aN79+4A9OjRgx49evDXX3/RuXNnQkNDefPNN+3Hv/XWW/Zx//77b3bv3s2QIUMAW6KAUqnk4MGDvP/++4AtRm/69Om89NJLDBo0qF7yNwViudOLWLNmDVu2bMFoNPLOO+/Qs2dP4uLiCA8PJyYmhoULF2KxWPjhhx8cgvcjIiLIzMzEaDQCoFAouPrqq3nppZfsmS2ZmZn2df9LLrmEn3/+mcOHD1NaWurSU1gFzz33nENWZdV/tS1ZTJw4kT/++IMtW7ZQUlLCO++8w+jRox1c2hWYTCbuvvtutFotr776qt1jVcHGjRs5e/YssiyTnp7OG2+8wahRo+z7Z82aVeeyxldffUVGRgb5+fn2IOYKxo0bx/r16/n6668djBu9Xo9GoyEsLIzS0lLeeustp+NXxFDs27cPg8HgEMdR19+pMdeuXbt2dO7cmffeew+DwcDy5cs5cOAAY8eOBWDZsmUcPXoUq9VKbm4uL7/8Ml26dCE0NBSwxXd8//33lJWVUVZWxrfffmuP7dixY4f9e1pWVsbcuXM5d+6cXeEJvI/k5GTGjRvnkAk4YsQIjh8/zqJFizCbzfz2228cPnyYCy+8ELB5/RuTJFQfnN1/td03YNN5p0+ftr/v0aMHgYGBzJ07l7KyMiwWCwcPHmTnzp2ALRFo7ty5FBQUkJGR0aCMyMWLFzvVebU9QE6aNIn58+eTmZlJZmYm8+fPd5qFPnr0aA4dOsTSpUsxGAy89957pKWl2Y2q77//3q4rDh8+zNy5cx0MhhtvvLHOWLHZs2dTWlrKoUOH+Omnnxx03sSJE1m4cCFLly51iC3W6/X4+/sTHBxMZmYmH330kdPxO3fuzJo1a8jPzyc7O5tPP/3Uvq+uv1Njrl337t3ZsmWLfUl97969bN261a63Tp48SV5eHhaLhTVr1vDtt98yffp0AO655x6WLl3KggULWLBgASNHjuTqq6+2x1AePHiQW265hSeffJKRI0fWen2bGmGkeRETJkzgvffes2cRvv766/Z9zz//PB9//DEDBgzg8OHD9O7d275v4MCBpKSkMHToUAYMGADAgw8+SHJyMtdccw19+vRh6tSp9qzIESNG8O9//5t///vfjB49moEDB3rsM1a45x944AEGDx6MXq93KP3w1FNP2YOUt2/fzh9//MH69evp37+/vRZaRXLE3r17ufbaa+nVqxdTpkyhY8eOPP744/ax0tPT6dOnT63yTJgwgZtvvpmLL76YpKQk+00LNk9DRUJGVUV2+eWXEx8fz7Bhwxg/fjy9evVyOn67du2YMWMGU6dOZcyYMdVit2r7OzXm2oHtSXH37t3079+fN954g3fffZfw8HDAZgzecsst9OnTh4kTJ6JQKByM9ZdeeokzZ84wYsQIhg8fzqlTp+zxfkajkeeee44BAwYwfPhw1q5dy9y5c+vlDRQ0HzNmzHBYbgwLC2POnDnMnz+fAQMG8NFHHzFnzhz7d+Rf//oXS5cupX///rzwwgtukcnZ/VfXfXPVVVdx+PBh+vXrxx133IFSqeT9999n//79jBo1ioEDBzrUsrrzzjuJj49n1KhR3HzzzUyaNMktn6cmpkyZwkUXXcTEiROZOHEiI0aMcAj2Hz9+PL/88gsA4eHhzJo1i7fffpv+/fuzc+dOh4fAbdu2MXHiRHr16sW0adMYPnw4999/v31/fXTeBRdcwOjRo5k6dSo333wzQ4cOte8bOXIkx48fJzIy0sErfuedd7J371769evHtGnTas14nTRpEp06dWLkyJHcfPPNDrqzrr9TY67dBRdcwF133cXdd99N7969ueuuu7jtttvsn2/37t1MnDiRPn368NZbb/HGG2/YvXKBgYH28JmoqCj8/PzQ6XT2h9b58+eTm5vL448/bv8dGj9+fK3XuamQZFfXhwRu4ZFHHiEmJsahJo6gcUyaNIlPPvmEsLCwGvePHDmSF154gcGDB3tYMoFAIO6/piUjI4N77rmHb7/9tsb9FWV19uzZg0olIp18BfGXErRYFi5c2NwiCAQCgUeIjY11aqAJfBex3CkQCAQCgUDghYjlToFAIBAIBAIvRHjSBAKBQCAQCLwQYaQJBAKBQCAQeCEtMnFgwIAB9aqiLxAIWgZnzpxh48aNzS1GkyD0l0DQ+nCmw1qkkZaQkMBPP/3U3GIIBAIPUbW1lq8j9JdA0PpwpsPEcqdAIBAIBAKBFyKMNIFAIBAIBAIvRBhpAoFAIBAIBF6IMNIEAoFAIBAIvBBhpAkEAoFAIBB4IcJIEwgEAoFAIPBChJEmEAgEAoFA4IUII00gEAgEAoHACxFGmkAgEAgEAoEXIow0gUAgEAgEAi+kRbaFEggEjSMnJ4e77r4LvV5PZGQks9+bjVarbW6xBAKBoFUhPGkCgaAaJ06c4PSp0+Sacjl08BA5OTnNLZJAIBC0OoSRJhAIqlFSUgKAnCADUFpa2pziCAQCQatEGGkCgaAaFUYauvPeCwQCgcBjCCNNIBBUw+5J08kO7wUCgUDgOYSRJhAIqlFcXGx74X/ee4FAIBB4DGGkCQSCauj1ept28KvyXiAQCAQeRRhpAoGgGnq9HoVGAWrbe+FJEwgEAs8jjDSBQFCN4uJim4GmBBTCSBMIBILmQBhpAoGgGkVFRVjVVpBAoVFQVFTU3CIJBAJBq0MYaQKBoBoFhQXIaltmJxooLCxsXoEEAoGgFSKMNIFAUI2CwgJkjc1Is6qswpMmEAgEzYAw0gQCQTWKCotAY3sta2QKCguaVyCBQCBohQgjTSAQOGCxWNAX6x2MtPyC/GaVSSAQCFojwkgTCAQOFBcXI8syaMs3aMs9awKBQCDwKMJIEwgEDhQUlC9tlnvS0EJZaRkGg6HZZBIIBILWiDDSBAKBAxVGmqytzO4EkeEpEAgEnkYYaQKBwIH8/Hzbi/Llzgpjzb5dIBAIBB5BGGkCgcCBvLw824vyvp0V/9u3CwQCgcAjNKuRtnbtWsaOHcvo0aOZO3dutf1Hjhzh2muvpVu3bnz88cfNIKFA0PqwG2NVEgdAeNIEAoHA06iaa2KLxcJzzz3H/PnziYmJ4aqrrmLkyJGkpKTYjwkNDeXxxx9n5cqVzSWmQNDqyMvLQ6FRYFFYbBvKjbTc3NzmE0ogEAhaIc3mSdu5cyfJyckkJSWh0WgYP358NWMsIiKCHj16oFI1my0pELQ6cnNzK5c6AdQgKSWx3CkQCAQeptmMtMzMTGJjY+3vY2JiyMzMbC5xBAJBOTm5OVg0lsoNEkg6SXjSBAKBwMM0m5Emy3K1bZIkNYMkAoGgKtnnspF1jvenVWMlJyenmSQSCASC1kmzGWmxsbFkZGTY32dmZhIdHd1c4ggEgnLycvMclzsB2U8m+1x28wgkEAgErZRmM9K6d+/O8ePHOXXqFEajkcWLFzNy5MjmEkcgEAAlJSUYygzVjTSdLDxpAoFA4GGaLSJfpVLx1FNPccstt2CxWLjyyitJTU3l66+/BuC6664jOzubK6+8kuLiYhQKBZ9++im//fYbgYGBzSW2QNCisRtiuvN2+EFxUTFGoxGNRlPtPIFAIBA0Pc2aNjlixAhGjBjhsO26666zv46KimLt2rWeFksgaLVkZ9uWNM+PSasw2nJycoiLi/OwVAKBQNA6ER0HBAKBnXPnztlenOdJqzDa7PsFAoFA4HaEkSYQCOxUeNLOj0mrMNrs+wUCgUDgdoSRJhAI7Jw7dw5JI4H6vB26yv0CgUAg8AzCSBMIBHays7ORdDXUK1SDpJLIysryvFACgUDQShFGmkAgsJOZmYnFz1J9R3nXAbHcKRAIBJ5DGGkCgcBOZnZm9czOciw6C5lZonWbQCAQeAphpAkEAgDMZjMFeQXgX/N+WSeTkZlR806BQCAQNDnCSBMIBIAtHk2WZadGGv6Ql5OH2Wz2qFwCgUDQWmnWYrYCQWORZRmj0VhtuyRJojK+i1QkBcj+NS934m+73jk5OcTExHhQMoFAIGidCCNN4NO88cYbLFq0qMZ9Dz30EBMmTPCwRL5LZmZ5vNn5LaHKqYhVy8rKEkaaQCAQeABhpAl8FoPBwIrly+gUaqJHpMlh3+qzOpYs+U0YaS5gL6/hbLkz4LzjBAKBQOBWhJEm8Fm2bNlCaZmBiZ3K6BnpGCdltkr8vHsPOTk5RERENJOEvkVmZiYKrQKLqoYSHGA33jIyRPKAQCAQeAKROCDwWZYtW0aQBrqGVw9kHxBjRJZlVq5c2QyS+SaZmZnO49EAVKDQKoQnTSAQCDyEMNIEPklRURHr1v3JoJgyVDV8ixMCrLQPsbLkt8W2jEVBnaRnpGPVWWs9RtbJwpMmEAgEHkIsd3oAg8HAfffP5MiRw3Uee83VV/Of//zHA1L5NkuWLMFkMjMivnpmZwUj4kqZv/8Yu3fvpnv37h6UzveQZdnmSUuo3aC16qyiVppAIBB4COFJ8wBffvklu3ftpDAgicLgdk7/FatC+fSzz9i/f39zi+zVWK1WfvrxBzqGWkgOchI/BQyJNeKvlvjxxx89KJ1vUlxcTFlpmfOkgXLkAOFJEwgEAk8hPGlu5tSpU3zx5ZeYIzpgbD+s9oPNRgJ3/cAbb7zJnDnvo1KJP09NrF27lrPpGdzVvbTW4/xUNm/a0tWrOXv2LPHx8R6S0Peos0ZaBf5QWlJKcXExgYGBHpBMIBAIWi/Ck+ZGzGYzzz73HBYUGNtcUPcJKg2lbQZy8OABPv/8c/cL6IPIsswXn39GbIBM/2hTncdf2qYMCStff/21B6TzXezesYA6Diz3tNlrqgkEAoHAbQgjzY3MmzePgwcOUJo8BFlT16+fDUtEe8wRHfj000/ZuXOnmyX0PdatW8fBQ4eZ2KYEhVT38eF+MsPjDPy2+FexTFcLdqOrruXOck+bMNIEAoHA/QgjzU2sXr2aL774AlNUGpaI9i6da2g7BKs2iMefeEKUO6iCxWLhw7kfEBcgMzTOecLA+VzerhSsFubNm+dG6XybzMxMJKUE2joOLDfixPdSIBAI3I8w0tzAgQMHeOHFF7EGxWBsO9j1AVQaSlIuprBYz6OPPkZJSUnTC+mDLFmyhOMnTnJlez1KF765EX4yoxPLWLZ0KYcP151h2xrJzMxE8pegLu+kH6AQBW0FAoHAEwgjrYk5deoUM2c+gFHSUJYyChTKBo0j+4dR2v4iDh0+xONPPIHJVHf8VUtGr9fz4QdzSA21MKAesWjnM6ldGQFqmDXrXVE3rQYyMjOw+DnPlLUjgSJAFLQVCAQCTyCMtCYkIyOD++6/n6IyIyUdL0HW1BHgUweWsDYY2g1j65YtPP/885jN1SvrtxY++eQT8goK+b9UPVI9YtHOJ0Atc0U7Pdu372D16tVNLp+vk5GRgRxQP+PV4mcRMWkCgUDgAYSR1kSkp6dz5113kZ2TT0nqGGRdSJOMa47qiKHNAFavXs0zzzzbKg21I0eO8P3333NRgoEOIfXw9jhhZIKBtsFWZr37jlhCroLZbCYvNw909Tte9pc5m37WvUIJBAKBQBhpTcGpU6eYceddZOfmU5J2CdbAqCYd3xzXHUObAaxdu4Ynn3wSg8HQpON7MxaLhTdef50AlZVrOtReF60ulAqYmlZMTk4uH374YRNJ6PtkZ2fbloDrl4AM/pCXm9cqHxgEAoHAkwgjrZHs37+f26dPJ6egiJK0S5vcQKvAHNcdQ/Ig1q9fz8yZMykuLnbLPN7GggUL2LN3Lzek6gnSND6WLCXEwsWJZfz004/s3bu3CST0feyFbHX1vL46W9eHnJwcN0olEAgEAmGkNYKNGzdy1913U2SU0XeegDUg0q3zmWO7UtbhInbu3s2MO+9s8cHbGRkZzP1gDt0jzAyJrX/Jjbq4OqWUMD945eWXMBqbblxfxf49qmcIZUWttJb+/RMIBILmRhhpDWTBggU8/PDDlCkDKOk8AdmvaWLQ6sIS2YHSjmM4fvI0t067jYMHD3pkXk8jyzKvvfYqVrORmzs1LFnAGf4quDmtiOMnTorODrhupIlaaQKBQOAZhJHmIhaLhf/973+89dZbmIITKOk8vt7dBJoKa0gi+s4TyC8xMGPGnfz5558end8TLF68mC1btnJdSjFROmuTj98r0szQOANffPF5izV060tWVhYKjaL+nXx1lecJBAKBwH0II80FiouLeejhh/nuu+8wxXShrONoUGqaRRbZPxx958soU9s6E3zxxRctpv5XZmYm/5s1i85hZkYmuG858v86lhKksvLySy+16jp02dnZdTdWr4oaJLVEdna2+4QSCAQCgTDS6supU6e47fbb2bx5C4Z2Q22dBKTmvXyyxp+SzuMxh7dj7ty5vPDCCz6f+SnLMq+//joWYxm3dtbXqz9nQwlUy9zUqZgjR4/yxRdfuG8iLycjMwOrnwveSgkknSQ8aQKBQOBmhJFWDzZv3syt06ZxOj2L0k6XYI7u1NwiVaJQYehwEcbEvixfvpwZd97p0x6OpUuXsmnTJq7poCfav+mXOc+nb5SJwbEGPv/sM44ePer2+byRrKys+md2lmPxs5CVLYw0gUAgcCfCSKsFWZb58ccfefDBBymRtei7TMIaHN/cYlVHkjAl9KYsdTSHjhzllltuZd++fc0tlcvk5uYy69136BhqYXSS5zyCN3YsRae08MrLL2GxNLxYri9iMpkoLCisdyHbCmSdLIw0gUAgcDPCSHOC2WzmzTff5J133sEYkoi+ywRkv6DmFqtWLOHJ6DtPJK/UxJ133sXKlSubWySX+N///kdpiZ5bOhe7dZnzfII0Mjd2LGb/gYMsWLDAcxN7AfZaZy4aaeggPze/1Rm1AoFA4EmEkVYDxcXFPPDAA/zyyy8Y43piSG2+BAFXkf3D0Xe5DINfOM8++yyffPKJTyQUbNmyhRUrVjAxuZT4APcvc57PoBgT3SPMfDj3A59eLnaVc+fOAS4Usq2gvKBtfn5+0wslEAgEAkAYadXIysrijhkz2LZ9B4b2wzG16U+TFunyBGodpZ0uxRSZwrx583jttde8uoWP2Wzmv2+/RYy/zMS2Zc0igyTB1DQ9JkMZc+bMaRYZmoMKI60hy50O5wsEAoGgyRFGWhVOnDjBbbffzolTZyhNG4s5qmNzi9RwFEqM7UdgjO/F4sWLeeyxx7w283PBggWcPHWaG1L1aJTNJ0eMv5VL25SyfPly9uzZ03yCeJDGLHc6nC8QCASCJkcYaeUcOXKEGXfeRW6hnpJO47GGJDS3SI1HkjAl9cPQdggbNmzk4YcfprS0cU3Kmxq9Xs8n8+fRLdxM78jmr1V2WdsyQv3g/dmzfWKZuLHk5OTYtICrq/l+tv+EJ00gEAjchzDSgGPHjnHX3fdQVGZC32k81oCI5hapSTHHdMbQfjjbtm/nwYce8iqP2g8//EBhUTHXpJR4xaqynwomJZewc9cutmzZ0tziuJ2cnBwUOgW4eu2FkSYQCARup9ZGMBs2bGDZsmWkp6ejUqlITk7m6quvJjk52VPyuZ3MzEzunzkTvdGCvtN4ZL/g5hbJLZijUpElBTv/Wc1zzz3Hc889h1LZjGuLQGlpKd99+w19Ik20D/aeLMELEwz8etKfTz/5hP79+ze3OG4lJycHq7YBiRoKUPgpyMvLa3qhWgitQX8KBAL34tST9sYbb7Bw4UJ69uyJWq0mMTGRNm3acM8997BkyRJPyug2SktLeeDBB8nNL6QkdUyLNdAqsER2wJA8gD///JPZs2c3tzgsXbqUomI9E9p61xKsWgFjE23etAMHDjS3OG7l3LlzyNqGLevKfrKISXNCa9CfAoHA/Tg10tasWcPLL7/MpEmTeOutt9i+fTvXXHMNn376Ke+9954nZXQbs2bN4sTxE5R2GNniljidYY7thimmC99//z1//fVXs8khyzI//fgD7YOtpIZ4jxetggsTDGhVUouvm5aTl4Ps1zAjzaqxkpMrjLSaaA36UyAQuB+nRpokSfYaSFlZWVittiWRkJCQFhFQvX79en799VeMcd2xhCY2tzgexdhmAHJABC++9FKz1bk6evQox0+cZHh8qVfEop2Pvwr6R5Wx+o9VGI3ua/LenFgsFgrzC+3xZa4i+8mcyxExaTXR0vWnQCDwDE6NtNtvv53Jkydz8803c/3113PHHXcAttY9nTp5Ue/KBmCxWHhv9mzwD8OU2K+5xfE8CiWl7UdQVFTE119/3Swi/PHHHygkuCC6+TM6nTEoxoi+pJTNmzc3tyhuobCw0GYwNNBIww/y8/KF0VEDLVl/CgQCz+E0cWDcuHEMHjyYU6dOkZycTHCwLV4rPDycN998s0kmX7t2LS+++CJWq5Wrr76aadOmOeyXZZkXX3yRNWvW4OfnxyuvvELXrl0bPe/KlSs5feoUZamjQNE6E1xl/3BMER344YcfmTJlCmFhYR6df+uWLXQIsRCs8d4f+C7hZtRKiW3btjFkyJDmFqfJqQj6b2hMGlowGU2Ulpbi7+/fhJL5Pp7QnwKBoOVTa3ZnaGgooaGh1bYfOXKEDh06NGpii8XCc889x/z584mJieGqq65i5MiRpKSk2I9Zu3Ytx48fZ9myZfzzzz8888wzfP/99w2ec8SIEaxduxYAtVpDcJ5ESEpvQlP6AJC7fyNFx3eji04issdFKFRq9OlHyNm9DpUukKjeo1AHhGIszCF7x0qsJgMR3YbjH5OM1WQge8cqynLPEtLBN8bM2fMXWM1ERETw888/M2nSJJev6VdffcXKlSuJioqioKCA0NBQZsyYQWJiIunp6fzvf/8jLy+PqVOncsEFF2AwGJg1axYLFi5kVKo/YPtx/31fIasOFpEapeXG/uH4qRXsOlvKd9vzCNQo+NcFEcSFqMksMvHpxlwKyixc3SuUXon+GMxWvtycy77MMkakBDGhW0iTjXnmaCY///QTd911FwA///wzixcvpkuXLkyfPh2dTse2bduYN28eSqUSlUqFXq/n+uuvZ/jw4Q3+rnoCe2amtvo+/Vk9WX9nIakkYofG4hdRg7vNr3IcYaRVx536UyAQtBLkBjBixIiGnObAtm3b5Jtvvtn+fs6cOfKcOXMcjnnyySflRYsW2d+PGTNGzszMrHPsyZMnV9s2bdo0GajxX5vRU+WE4dc6bAvrNFDueP2TMgqlfZsmJEruPn2WrA4ItW+TlGq507+fl4Pb9WzEmO8225hV/7nKiy++WOM48fHxcmZmppycnGzfplKp5A0bNsjXXHONw7FPjI2VXxgf57Dt4rQg+bfbO8gapWTfFhOkkrc+mCYnhqorx1QgL7i1vTyha7Dbx3zjjTfk2bNnO2ybOHGivH37dlmj0VS7BpIkyUuXLnX5mnqSFStWyMOGDZMHvzBYHvxu5b/ej/eWFWpF5TUJVMn9X+7vcMzgdwfLgx8dLA8bNkzevXt3c3+UGu95b6Uu/elLn0UgEDQNzu57p560F154wZlRR2FhobPT6k1mZiaxsbH29zExMezcubPWY2JjY8nMzCQ6Otrl+T788EOn+87tXI3F6FgGIu/ARlAqwVqZeWgsyCZz82+Y9Pn2bbLFRPb2lRQe+6cRYy5p1JiSUtXgMasydepUPvnkE6f7z2fu3Lk1bj979iyvv/46J06csG8zm8188MEH1TyhX2/Nw0/tmDmw4kARUYEqjJbKZbjMIjMfrD/H6fzKGDazFb7cksvivY7fR2djRtc05l851cb8anP1MT/44AMCAwMdti1atIi4uLgaEwtkWebjjz9mzJgx1fZ5C/akkfM8adlbsrGaKmunmYvN5O7MJWZQjOOB5eeJWmnVcbf+FAgErQOnRtqPP/7II488gkZTvV/Mr7/+2uiJ5RqCjaXz0vzqc0x9USgUWCw1l3pQ+gVUa6KuUGlR+QVWO1YVUL2Wmso/CEmpQrZUNjH37JgBDR6zKt26datxuzPCwsIcDLGqxMTEVNsWERGBv78/er3evi3YT4Gf2jEuUKuSCPevXmg3MrD61zVMp0SnVlBirDQqnI0ZVtOYAdW3hdYwpr+/f7WlK61WS2RkZLXz7bJ5OM7PVfLz822dBs67xVW66tdZ5V+Dqig30goKCppcNl/H3fpTIBC0DpxGzXfv3p3U1FQmT55c7V9AQHWjwFViY2PJyMiwv6/JQ3b+MRkZGQ3yooGtcGpNKFQa4gZNIn7IlUhKtX173ODLiek7FnVQuH1bSIfexPS9hKA2XezbtKHRRPcZQ8wFE5ptzOgGjnk+DzzwgNN9NfH888+jVqurbb/44ou57777mDChcq6kpCTuvfdenn76afs2jVJi5sgY7rsw2sHzNWNYFDcPjCQprHLsoe0DuGVgBKM6Btm3xYeo+c+gSO4ZEVX5OVVNMOZgxzElSeL222/nmWeeQaer7ET+2GOPcc8999CuXbtq1yAqKoqZM2c6uXLeQX5+Pgpt9ZZQ0QOj0cVUfs6g9kGEdavB4NRWjiNwxN36UyAQtA4kuSZ3FTbFq9VqHX6UmhKz2czYsWP55JNP7IkDb775JqmpqfZjVq9ezRdffMGHH37IP//8wwsvvMAPP/xQ59hXXHEFP/30U7XtZ86coVOnTuj1ejr0Hok5Po2gNp1R+9u8TqbifIpO70cX1QZdRDwAVpOBwuO7UPoFEpiYhiRJyLKVopP7sJoMBLfthkJle1ouyT5JWU66z4xZfPYwOduXExUVRWZmput/xPJr+ueff9K5c2eOHz9OUFAQF110Ufn8MmvXriU3N5exY8fag8tff/11Pv74Yz68FJLLjaZzxWb+OlZMx2g/OsXYItLLTFbWHikmQKNgcLsA+5gbT5SQX2phRIdAdBrbc8bBrDL2Z5YxoG0AMUFNN+ayIybW5sXywQcf0LNnT7Kysli1ahVdu3ale/futjHLyli6dClBQUFotVrS09MZO3YsQUGVxp838tRTT7F2+1pMY6qXQbGareTvy0ehVhDSMQRJUYMHWwbVzyquueoaZsyY4QGJnePsnm8uGqM/ve2zCAQC9+Psvne63FlSUlJjZlJToVKpeOqpp7jllluwWCxceeWVpKam2ut2XXfddYwYMYI1a9YwevRodDodL730UqPmTEhIoKioiAceeIDN/+ylOK0/SJXORHVgKOGdBjqco1BrCU11rKUmSQqCk6uXAvGPaoN/VBuHbd48ZmRiezoHlHL//fdXO6a+JCQkMGXKFAB69ux53vwSI0aMqHZOWloa0dHRBOryscWm25YyL+se6nCcn1rBmE6Oy7aSJDGwbXVPRMdoPzpGO2YgNsWY+XIQG3Zp7PFo0dHR9s9rH9PPr0GZsc1Nfn4+FnXNIQAKlYLw7uE17rMjgeQniRirGnC3/hQIBK0Dp0bajBkz+PnnnwG46667mDVrVpNPPmLEiGo/4tddd539tSRJDstjTcXEiRPZtGkTyrxTWMJbb7NjVdY+NFotF198sUfnrYjVyjcoiPDzvpZQVck32oz4kJCQZpak6cnNy60Wj+YqskYWMWk14An92ZwsW7aM9evXExgYyN13341WW0MdF4FA0GicxqRVXQU9deqUR4TxFEOGDCE+PgG/05scsiJbEwr9OdTnDjFxwoRqWYvuJikpCYD0kupB+95Gul6Bv05HRETL6+1aUFjQ8EK25VjVVvIL8ptGoBZES9afAJ9/8QV/rF7DokWL2LdvX3OLIxC0WGrt3VnT65aASqXi/vvvg9IC1Gd3NLc4nsdqxe/4ekJCQrj55ps9Pn1CQgJqlZKTxd5vpJ0sVtG2bdsWdw/IskxxYXHjPWlambx8UYLjfFqy/gTIycnBEhhlfy0QCNyD0+XO/fv306dPH2RZxmAw0KePrdq9LMtIkq1Vji9zwQUXMHr0aJavWIElKAZrSOtpsq45uRGpOJv7n322WYLb1Wo1nTp14uCpnUBpncc3FyYrHC1UccXYHs0tSpNTWlpqK0nT2FUqDRRmipi082nJ+tNgMFBcVIQlri3Kokyys7ObWySfZ9euXcyfPx+LxUJSUhIzZ85skca9wHWcGmmtwYX9wAMPcPDQIU4eWY2+y0Rkv5YXd3Q+quyDqDP3cPXVV3PRRRc1mxw9e/Xm6z17KDFDTSW4vIFD+SpMVujRo+UZafY4skZ60tBAib4Ei8WCUun9nlFP0ZL1Z1ZWFgCyLhxJpWlwZrigkuXLl7Nj2xYi/Sxs376dm2++mfDwOhJ3BK0Cp8udf//9t/31+TEVy5Ytc59EHkSn0/HSiy8SqNPgf+B3JKO+7pN8GGXeCbTH/qR3795Mnz69WWUZNGgQFhl2nqteZ81b2JatRq1W0bdv3+YWpcmpyMiUG9vgXlu+dFpc3ARStRxasv5MT08HwOoXhKwN4uzZs80ske9z8uQJkoMs/CutBMBpkXBB68Opkfbaa6/ZX999990O+95//333SeRhkpKSeOvNN/GTzPgf+B1MJc0tkltQFJxBd3gVaWlpvPzyy6hUzeu+6tKlC2EhwWzObqwrxz1YZdhyzo++ffu2yObhdk9aEyx3AqIMx3m0ZP1ZYXTK2mDMmkBOtsDECE8iyzJHDh8mMcBMUqAtke3o0aPNLJXAW6hXduf59W6d1L/1WdLS0njt1VfRWEoI2PcbkqFledSUeSfxP7iMtsnJvPH6615hdCiVSi4cOYpt57SU1Nylqlk5VKDkXCmMGuXZ8iSeoqioyPaiCUpwgGgNdT4tWX+eOnUKSaVBVuuw6kJJT0/HbPbCm9hHyMrKoqCwiORAC6EamRCtxIEDB5pbLIGX0KDszpYY0NirVy/efust/DDgv38xUlnL8Awoc47id2gFKSkpzJr1rlfV+xozZgwmi8ymTO/zpq1L1+Kn1TBs2LDmFsUtNGVMGlQx+gRAy9afR48exeIXCpKEVReK1WLh9OnTzS2Wz1IRv9g+xIwkQfsgA3v37G5mqQTegtM1r1OnTnH77bdXew202Buye/fuvPPf/zLzgQeR9v1KSccxWAOcN9D2dlQZe9Ce2EC3bt147bVXPV4PrS66dOlCcpskVp89zoUJxuYWx06ZGf7O9OPCi0d6hdfRHTSVJ00sd9ZMS9Wfsixz6PARLDpbOzpZZwtuP3z4MG3btm1GyXyXf/75B40S2gbZljpTQ81sP3ya/Px80bVC4NxImz17tv31+bW0mqO2lqfo3Lkzc96fzb333Q/7F1OSMsr3ynPIMurTW9Cc/YchQ4bwzDPPeGVFcEmSmHjZJP73v/9xskhJmyDvKCz8d6aGMrPMZZdd1tyiuI2CggIktVSLL72eaCvHE1TSUvVnRkYG+uIirJG2h1erLgwUSg4ePOjxziUthW1bt9AxxISq/F7sHGpbOt66dSujRo1qRskE3oBTI+2CCy7wpBxeRZs2bZj7wRzunzmT4weWYWg3DHNUat0negNWC5pj61CfO8SECRO4//77mz1JoDbGjh3L3A8+YOUZDTd1av6aabIMK07raN+uLV27Vu972lIoLCxE0jbBslt5cq5Y7nSkperPiqU5+wqDQoHVP5y9e/c2o1S+S1ZWFseOn+C6VJN9W/tgCwFqiU2bNgkjTeD8OXrFihV8+eWX9vdXX301o0aNYtSoUfz+++8eEa45iYyMZPZ779G7V0+0R9egPrPd9gvuzZiN+B1chvrcIf7zn//w4IMPerWBBraemCNHjWJdhs4rEgiOFCo5UaRg8hVX+nzsUG0UFhYiq5vg+yyBQqsQnrTzaKn6c8+ePUgKFVb/yjZplsBo9u3fj8lkquVMQU2sW7cOgF6RlddOqYAe4Qb+Wr9OJGQInBtpH330ESNHjrS/NxqN/PDDD3z++ed8/fXXHhGuuQkMDOSNN95g9OjRaE5vRXP8L5CtzS1WjUjGEvz3L0ZdnMGjjz7Kv//9b58xMiZPnozBLLM+vfmXZFec0uKv82P06NHNLYpbyS/Ix6ppou+yRnjSzqel6s9t27djDowCReVPhyUoFpPRKDISG8CaNauJC5BJCHC8F/tFGykoLGLnzp3NJJnAW3BqpJlMJuLi4uzv+/btS1hYGPHx8ZSWNv+ylKdQq9U8/vjjXH/99aiz9qE9tBKs3vV0I5UW4L9vEX5mPa++8gqXXnppc4vkEp07dyatYyorzuia1VlZaJTYkKXlkkvHtdiEgQoKCgoaX8i2HKvaKjxp59ES9WdhYSFHjxzBEhTnsL3ivS+3umoOsrOz2bHjHwbFlFXb1yvShJ9KYsWKFc0gmcCbcGqknZ+t9dRTT9lf5+bmuk8iL0ShUHD77bdz9913o8o7ge7AMrB4RzaiQn+OgP2/EqxRMGvWuwwYMKC5RWoQl0++gjPFEgfym295du1ZDWYrTJo0qdlk8BSFhYWNz+wsR9bI5BfkN81gLYSWqD+3bNmCLMtYQhIcd6j9kAMi2bRpU/MI5qMsXboUWZYZHFv9t0SrhP5RZaxaudJnjXpB0+DUSOvRowffffddte3ffPNNi+xlWB+uuuoqnnjiCVTFGfjvXwLm6k9AnkRRlIn//t+ICA7k/fdn06lTp2aVpzGMHDkSf52OVWeap2aaLMMf6Tq6d+9Gu3btmkUGT2GxWCjRlzS+20A5slYWnrTzaIn68++//0ZSa7EGRlXbZwpJYPfu3WLZu55YLBYW/bKQTmFmYv1rDjsYEW+kpLSUVatWeVg6gTfh1G3x2GOPMWPGDBYtWmTPctuzZw9Go5H33nvPYwJ6G2PGjCEgIIAnnnwS//2/U9LpElD5eVwORVEm/gd/Jy4mmnf++19iYmI8LkNTotPpGDN2LIsXLURvKiWgKYLaXeBAvopMvcR/JrbcshsVFBUV2areN5U9rIHCTFEnrSotTX+azWbW//U3puBEkKo/21tCk7Ge/YcNGza0+HjOpmDjxo2kZ2RyZXfnD/ppoWYSg2R++vEHxo0b5zMxxoKmxaknLSIigm+++YY77riDhIQEEhISuOOOO/j222+JjPTdAq9NwZAhQ3j5pZdQGfLx3/87mD279KkoyrIbaP+bNcvnDbQKxo0bh8kiszHT803X/0zXoPPTMmLECI/P7WmarG9nBVowlBkwGAxNNKDv09L05z///ENxUSHm8LY17rcGRiFpA1izZo1nBfNRvv32G8L8oF+U84xYSYIxCSUcOnyEHTt2eE44gVfh1JOWn58P2IK6O3fuXG17a6+EPHDgQF568UUefewxdIeWU5o2FhTuj6eSSvPwP7SMmKgo/jdrlk8qfGekpaXRtk0S6zOOMTLRc4av0QKbsm0dBnQ6ncfmbS4q4qWaKnGgateBqKjqS2GtkZamP1etWoWkVGMJSar5AEnCGJrM339vQK/XExAQ4FkBfYh9+/axffsOrkstsRewdcbQOCM/HAvgyy+/pHfv3p4RUOBVOLUqrrjiCiRJQpZlsrOziY6OBmxtQSRJYuXKlR4T0lsZNGgQjz/2GM8//zzaI6sxpIyyPf64CcmoJ+DAUoIDdPz37bdalIEGtg4EIy8ezfx588gzSIRpPbPkuStXTalJbjWFIysMhaaMSQObh04YaTZakv40mUys+uMPTKFtQOn8QdQc0QFT5l7Wrl3rcxnmnuSzzz4lQA0XJdTtedYoYWxiKd9v2sT+/ft9Ou5Y0DCc3nFVgxUvv/xyFixY4Al5fI7Ro0eTk5PD7NmzsZ7dgSnBTU87Vgu6w6vQYObNN/5LfHy8e+ZpZi688ELmzZvH1mw1F3vIm7Y5S01QYAB9+vTxyHzNjTuWOx3GFbQo/bl+/Xr0xcWYE4bUepw1MBp0wSxZ8rsw0pxw4MAB1q//iyval+Jfz4WX0Ull/HbKn/nz5/Hqq6+5V0CB11Gvzn0iYLF2rr322vKCt9tQ5p9yyxyaExuQijJ57LFH6dixo1vm8AaSk5OJj4tlxznPxKVZZdiZq2XgoMFe352hqbB70poqcUB73rgCB3xdfy5e/BuSNqB66Y3zkSSMEans2LGds2fPekY4H2Pexx8ToJG4pE39KwP4q2BcUgl//71BtN9qhTS2vbIAmxJ+6KGHaNuuLbpjf4KpaUtzKPNPoc7ax7XXXstFF13UpGN7G5IkMXDQYPbmaTB5oLnDiSIlhQZ8tr5cQ8jPz0dSSbX40V1EGGktloyMDDZt2oghsmONWZ3nY47qCJLEokWLPCCdb7Fz507+3rCB8W309faiVTA2qYxgLXwwZ44tM1vQanD6VZk/f779dU5OjsN7gJtuusl9UvkgWq2Wp558kltvvRXtib8xpDSRMWU2oDu+jjbJbbn11lubZkwvp3fv3vz0008cK1TSMdTi1rn259lugday1AnlRppfE3p3NIAkjLSqtBT9uWjRImTAHJVWr+NlTQDm0Db8smgRU6dORatt/lZv3oAsy7w/ezahfjA2yfUsaD8VXJ5cwmc7drBx40YGDhzoBikF3ojTRyO9Xm//d8011zi81+v1npTRZ0hJSeHGG29ElXMERVFGk4ypPrsD2VjC448/hkbTPIVePU3Pnj0BPNJ94GC+ivjYmBaXhFEb+flN2LcT7E3WhZFWSUvQnwaDgQULF2IObYOsDaz3eaaYLhQVFvLHH3+4UTrfYu3atezZu5cr2+nRKhs2xshEAzH+MnPen43F4t6HV4H34PRX8M4773R6UklJiVuEaQlcf/31/LJoETknN1HSZWKjsj0lQxHazL2MHTu2VWX1hIaGEhsTxYmiM4B7a28d02voOaibW+fwNnJyc+wZmU2FrJXJy8tr0jF9mZagP1esWEFRYSGmTkNdOs8aHA/+YXz3/feMHTvW52PyGovJZGLO+7NJCJQZHtfwZCiVAq7poGfWruP8/vvvjB8/vgmlFHgrtQYZZGZmsmvXLoxG2xcrJyeHt956izFjxnhEOF/Ez8+P/9x8M1JxForC9EaNpU7fjUKCW265pYmk8x1SUtM4oXev51BvkjhXYvOAtiby8vKa3Eizaqw+25PSXfiy/pRlmW+++RY5IAJrcFzdJ1RFkjDEdOXwoUOiCCvwyy+/cOZsOtelFKNsZBT4BdEmUkIsfPThXNHTs5Xg9CvzySefMGnSJF544QWuvfZafv75Z8aNG0dZWRk//fSTJ2X0OcaMGUNwcAiajF0NH8RsRJNzkJEjR9prLLUm2rZtS5YeLG5MHsgosX39k5OT3TeJlyHLMgX5BU1XfqNiXD+ZnNycph3Uh/F1/blx40ZOnDiOMaZbg1YDzJEpSBodX331tRuk8x2KioqYP+9juoab6RlhbvR4kgTXp+rJyc3jm2++aQIJBd6O0+XO7777jt9//53Q0FDOnj3LmDFj+OKLL+jVq5cHxfNNtFotEydO4MuvvgJTKahdr2KvyjuObDZxxRVXuEFC7ychIQGLDOfKFMQ4aUDcWDJLbUZaYmKiW8b3RoqLi23xLE3dblYL+Wfzm3hQ38XX9eeXX32FpA3AHNG+YQMoVBiiOrNx4waOHj1K+/YNHMfH+fLLLykqKub6ASVNVue8Y6iF/tFGvv7qKy677DIiIiKaZmCBV+LUk6bVau2tS+Lj42nbtq3PKBhvYOTIkSDLqPJONOh8Ve4xoqKj6dKlSxNL5htU9CPNNbivSkxumW3s1uSptC9JNrWR5gelJaWif2c5vqw/9+7dyz87dlAW0w0UDYxyx5ZAICnVfPXVV00one+QmZnJ999/x5BYA8lBTRvof21KKSajoVrWsKDl4dSTlpGRwQsvvGB/n5OT4/D+iSeecK9kPk5KSgoRkVFkFpzBHO1i0L/Viqoog2GTJrbaoNvw8HAACgzu+/z5RgV+Wi3+/v5um8PbqDDSZL8mrrXkVzl+XJyLMUwtEF/Wn1988QWSWos5un5lN5yi9sMY1ZEVK1dyyy23EBsb2zQC+gjz589Htpi5qkPT1s0EiPW3MjKhjF9//ZVrr72WpCQnPVUFPo9TI+2hhx5yeN+1a1e3C9OSkCSJ3r16svLPvzHIsktxHYqSHGSLiR49erhRQu8mJCQEgGJz7dctr8TMF5tzyS2xMLlHKD0SdJzINfDk4nT2ZpQxpH0Az42LJ0RX3SNQbJIIDg5yi/zeit2T5oaYtIrxhZHmu/rz+PHjrFu3DmN8L1A2PnHHFNsdddY+vvnmG+69995Gj+crnDx5kt+XLGFMUhmRutrDNX7Zlc/GEyX0TNBxVc9QZOD1lZks2JlPTJCax8bEMqBt9Yb1l7crY226jnnz5vH000+76ZMImhunRtrkyZNr3G4wGBz60gmc06lTJ1asWAHmMpfi0hQlOfbzWysBATalVFqLkWa2yFw17xiHs21LbJ9tyuHbm9rz1G9n2ZNue3pdsLMAhSTx9hXV485KzRKBgfWv/9QSsBtprodJ1k4VT5rAd/XnV199haRUYYptGqNS1gZiikhh0a+/MnXqVPsScEvnk08+Qa2EiW1r96K9szqLt/7IAuCLzbAnvZQ2YRreX3cOgPRCM//56gQbZqYReF6BtRCtzJikEn5dtZJ//etftGvXzj0fRtCs1Cvgx2KxsGbNGh566CEuuugilixZ4m65WgQVLmhFmWuNp6WyAlQqlT0uqzVSUancYHFupG04rrcbaABmK3y28ZzdQKtg/dHiGs83WCT8/Jo6OMu7ycnJsd31Td0aVRhpTvEV/ZmZmcmy5csxRnZsULKTM0xxPTAZTfz4449NNqY3c/LkSVauXMnohFJCNLWHFXy+2fF++XprHuvO01dFBiv/nKm53Ma4Nga0SonPPvuscUILvJZaS7pv3ryZRYsWsWbNGnr06MG2bdtYuXIlOl1TP4a3TCqWfSSDnjJTBjm71iAplUT2uAhNsPOMHIVRT2RUNEplw4N2fR1JklAplVhq0XEB2urPGGH+KtqEqTmZZ7Jv6xZX8/fVIoNa3Tq6OFSQm5uLQqfAIjVxxfLy5dOcHFGGowJf05/ff/89VquMKa57tX1Wk5Fzu1ZjyM0gJKUPwW3rXwBa1oViDmvDDz/+yHXXXdfiY0C/+uorVAqZS+vRRD1QqyC7ik3mr1HQLU7HqoOVG9VKiY7RNT9MBmlkRiWUsmTVKm655RYSEhIaLb/Au3DqSRs+fDhvvvkmffr0YfHixcyaNQutVuu1CsYbqXDtG/Mz2P/F02Ru/o2MDYvY/8UzmEtr9u4ASKYywsPCPCSl96JUKrDKzj1pvRP9uTitMqYsMkDJjReE898rk0gOtxlfPRN0PDuu5hgpq0yrM4RzcnKwat1Q0kQBCp1CGGnl+Jr+LCoqYuEvv2AOb4+srR6neXThu5xe9SXZO1Zy+IfXyd33t0vjm+J6oi8u5rfffmsqkb2SnJwcli9byvC4MkLqUTD6/ouiUVRRcfeMiOL2IZGM7RyMJEG4v5LXJiUQFejcn3JJmzIkSeb7779vio8g8DKc/uXHjBnDypUrWbJkCUqlklGjRrXaTMOGUvHEmHNkF1Zj5VOVuaSQ/IObiexZcxN2yWoiIKBlP23WB6ssU9c37qPr2rD+qJ6/jhWzYFc+F793mBEpgSy4pR1qpYIgP+dGmIStuGtrIvtcdtNndpYj+8nCSCvH1/TnokWLMJSVYUqp7kUzFGRTeNyxMHf2jpWEdx5U7/GtQdFYg2L55ttvmTx5cot9OPr1118xmS1cUs8m6pd1D6VHvI61h4tZtKeAZ3/P4KO/c3h5Yjz/vSIRjVJCpaz9exOmlRkUY2DJb4u59dZb7fG8gpaBU0/aE088wapVq5g6dSobN25k7Nix5Obm8ttvv/lMg+DmRqWy2cAKVfUAIIXGeSyUAhm1uqmDhnwP2SojSbUbFJIk0T/Zn6+35nEm34wsw+pDxby2MqtWAw1AIYHV6saWBl7IuXPn3GakWTVWzuWcc8vYvoYv6U+z2cwPP/yIJTgOa0D1MAyFSguS40+FUu16LKcxtitZmZmsW7euwbJ6MxaLhV8WLqB7hJm4gPrrlbYRWg5mG9h0ogRZhtP5Ju784TQKiToNtApGJxooLTOwfPnyhoov8FJqTRyQJIlBgwbxwgsvsGrVKt58801WrlxpK9QqqJOKJ+eodp3QhlXWCNJFtcFiKCX/8DZkaw2xQbLs1U/dnkCWZcwWC6oaLoPFKvP+n9lc9fFRHv3lDP+cKSW3xPE6Ogu0rYpSAWazqc7jWgomk4niouKmz+wsR9bJnDsnjLQKfEV/rlu3jnPnsjHF1hxnpg4IJrrPaPt7hUpDYJvOnNu5GpO+/klRlrBk8Avixx+9vy1WQ9i6dSvZ53K4ML7uWLT1R4u56csTTP3iOOuOFLPzrKO+Kii1cCK3/s3Y2wdbSAqysmRJy15Obo3UmjhQFbVazciRIxk5ciRlZU1fnK8lYrHYDAeV1p/O/36egqP/YC4t4uyfP3JqxacABCZ1JvWahx2MMllS2M9trZjNtj53qhoeI2atyebt1ba09c0nS9iTXkpcsIr0wsreeINqqCt0PipJpsxYf0Xo67it20AFfpB/Kh+r1YpC4b5OEb6IN+vPBQsWgF8QllDnBVETL7qesLQLKMvLIHfv35xd+x0ASu23dJzyOLqoerRWkxQYozqxY8dmTpw40eJ65q5YsQJ/tUTvyNof/A5nG/j3FycwlWdF/XmkmMk9Qh0eLKODVLSPrH8xQ0mCYbFlfLVvP2fOnBEJBC2IOjXprFmzHN6//fbbfPbZZ+Tl5blNqJaCyVR+s0oKFCoNYR37Y8jNwFJWmTRQfGofxaf2O5wnSwqMrch4qImK9kJqRfWlud/2Oj69/3O2jBcnxNM7UUeoTsnVvUN5YFTd5UvUCjAYvOsH051UxIu5a7kTHVgtVgoKXCs505Lxdv159uxZtm3bZiu7IdX+cxAQn4IuMomiE7vt2yyGErK2Lq33fKYo2zwtLYHAbDaz7s+19IksQ1NHuN2KA4V2Aw1spYMSQtVc2yeMUJ2S3ok6PpzSBnU9lzoruCDa9nuzZs0al+UXeC91GmnnV8ru3r07KpWKl19+2W1CtRRKS21PRrKyMr7Maq5ufJ2/zapQUVJa93JdS6bSSKu+r02YY9mMIK2ClCgtD4yK4a/7OvLG5Yn4a+r25GiUcqvqNWkP6nfXcme58SeSByrxdv25dKnNwDJHptbreKulupeoJp3mFLUOc2gSvy1ZYveWtwR2795Nsb6EvlF1h08khVUv+9MmTMPkHiH8Mq09C27tQK9E1xPHInVW2gZb+fuvv1w+V+C91LnceX78xMUXX+w2YVoaJSUltheKSiMtssdF5OxZh1weC+UXEU9w8nnVvRVq9PoST4nplVQYT1plda/PQxfHsDejjDMFJvzUEpd0DmbkrEOYrRDsp+CT/2tL36S6lZxWAYaS1uOxtMeLuXG5s2KelJQUN03iW3iz/pRlmaXLlmMJjkPW1q/zRkBcB/xj21GSccy2QaEkqpdrMXbmyBQKDq1kx44d9OvXz1WxvZLNmzejkKBreN1G2thOwYzvGsziPYUADO8QyBsrMzlbaEaSYPqQSB4e3bA+p93DDSzevZuSkpIWX4+uteDUSMvNzbU3uQZYuHAhu3btIjU1lWuuuabVB7bXh4osLrlKDzz/mGQ6/d+z5O77G5UukIhuw5CUjn8GWaVBr2/dldsrlns1NRhpHaP9WHtPR/ZnlREdqGLU/2wGGkBhmZXXVmTw7U3t65xDrZQpK2tlnjQJ9xlpuirztHJ8QX8ePXqUs2dOY247pN7nSJJE6jUPk7P7T0zF+YR1GoB/tGuxZZbQJCSVmj/++KPFGGn//LODtsEW/OsR5a1SSsy+pg0nc41YZZn3151j7RFbCIwsw/vrz3FD/3ASQ10vtN05zMyi41b27NlD//79XT5f4H04XRP6z3/+Y389e/ZsfvnlF7p27cr69esb7arPz8/npptuYsyYMdx0001OY1geffRRBg0axIQJExo1X3NhN9LOK8Ghi0wgYdhVxPS7BJVfDQHuSg0lJd6Vpu9pKoKrna1aqpQS3eJ0aJQShWWO6e4ZhfVbRtEoZIxGY6uplZabm4vCT0GdxecaimgNZced+rOpWL9+PVCedekCSo2O6D5jSBh+jcsGGgAKFabgRP5ct65FlMCxWCwcPHCA1GDXMsXbhGtoG6Els8jxPFmGrKKGLQWnhNjO27dvX4POF3gfTo20qj9cy5cvZ9asWUyePJk333yTv/92rdr0+cydO5dBgwaxbNkyBg0axNy5c2s87oorruCjjz5q1FzNib0ektLxiagk+ySn13xLxsbFNXYekJVqykpLW3WGZ4UnrabEgaqE+qu4MNVxqWZyj9B6zaFR2ArmtqTYmNrIyclxX9IAgBIUGtF1ANyrP5uKjRs3IQdGImtcWxazGErI3LyE06u/Rl+x7OkiltAk8vPyOHasYed7E6dPn6bMYKRtUMP09eSeoQ7v20do6JnQsMBRfxXEBMChQ4cadH5L4Ny5cxw7dqzGf74Yg+zUOVtWVsbevXuxWq1YLBb7+rZarW50ev3KlSv5/PPPAbj88su58cYbefDBB6sd179/f06fPt2ouZqT4uJyF3YVI60k4xgHvn4B2WIzDHL2/Ennf7+AosqSZ8XxJSUlBAVVb9HSGqg00uo+9n9XJTFn/Tn2Z5YxvEMgN/YPr/skbMudYMvCbQ3Fg8/lnHNPS6iq+AlPGrhXfzYFZWVl7Nm7B1N017oProIsWzn43SuUZp4AIGvrMlKveZigpE4ujWMJjgdg27ZtdOjQwaVzvY2TJ08CkBDYMCNtUvdQlJLEL7vyiQ/RcPvQSJSKhru743VGTp443uDzfZktW7bwwIMPYLXUrOfatmvLh3M/RKutf3mT5sapkRYVFWV3y4eEhJCVlUV0dDR5eXmNbumRk5NDdHQ0ANHR0S1WqVckDlTN7jy3a43dQAMw5KZTdGIPIe17Vp5YfnxrNtIqypfU5UkDCPJT8mA9Sm6cj1pynKulk5OTgxzk3qVdi9YiCtriXv3ZFBw8eBCrxYI10LX7Rn/2iN1AA0C2cu6fVS4babI2EMkvkL1797p0njdy5swZAGJ0DX8AmtAthAndQppEnlh/K3vPnkVuZUXR8/LyeP7F55G1MtZuVuTzfjukAonj+44ze/Zs7rvvvmaS0nWcGmkVnq7zCQ4O5ssvv6xz4KlTp9aorO+99976S+fj2ItWVsnuVKiqW/AKteO2CqOutBWX4ahYgnSxVJBLKBWVnrSWjtVqJT8/H6LcO4/sJ4vWUDRef7qbiuUwa2CkS+cp1NWD2c/XX/XFpItg/4EDDTrXm8jOzkarkvBXeUdsa4SfFYPRRFFREcHBwc0tjkcwGAw89vhj5OfnYx5phtDqx8hJMlaLlZ9//pkOHTpw2WWXeVzOhlDvjgMVKJVKdLq618s/+eQTp/siIiLsT5ZZWVkOWVAtCfv6t6LyyTmq98Xk7vsLc4kt/ToouSuBiWmOJ0oqx/NbIRXxeMp6eNIaSsWKQksIXq6LoqIi2xKAu738fpB3yjsKtXoj9dWf7ubEiRNIKi2y2rV4NP/oZEJT+5F/aAsASq0/0X0vaZAMVl0Y6Wd3+ny4QX5+PiFaGW9xWgVrbPosPz+/VRhpRqORJ596kj2792AZaKnRQKtA7i5DEbz55psEBAQwatQoj8nZUFw20gAmT57Mzz//3OBJR44cyYIFC5g2bRoLFizwiQvVEIxGo81Aq3L3akOj6HLzKxQc2YFKF0hw2+7VXdLlRl1r7jrgiYzLiqveGrI77RXu3VV+owItGMoMlJWV4efn7sl8k8bqz6YgPT0dqzaIhlgW7S6bQdGJPZiK8wlp3wuVf8NCMmRtEFarlezsbOLj4xs0hjdQVFSEv9J7HvQqPHpFRUXNLIn7KSkp4cknn2Tz5s1Y+1jBeWczGwqwDLSgXKfkueeeo6ysjPHjx3tE1obSoAjWxiqYadOmsX79esaMGcP69euZNm0aAJmZmdx666324+6//36mTJnCsWPHGD58ON9//32j5vU0FosFqYZWKyq/ACK6DiGkfU+kGoKI5XLF2Ro8PM7wRCxFyzfNKsnPzwfc2BKqAj/H+QTVaW4DDSArKxuLi160CiRJQXDb7kR0G9ZgAw1A1tjKD/l6DGNZWRlahffoam35wk1LX4nJzMxkxp0z2LxlM9Z+VuQO9dRtKrAMtSBHy7z66qt8+OGHXv1bW29PWnFxMcePHycpKYmQkMYFOIaFhfHpp59W2x4TE8OHH35of//WW281ap7mRpblBtakkirPb6VUZMBZZccLuOmEnsd/PcuxHCNjOgXx6mUJBPk1LBC74vJ6Q7adu7F70ty83ClrZft8sbENq5reEmlK/dkUFBUXI6uaVw5ZZYtvs5cq8lEsFgtKqX66urDMwsMLz7D8QBHtIjS8OCGeC5JrqJXZCCpkacklnLZu3crTzzxNUUkRliEWiHNxgHJDTdom8fnnn3Po0CGeeOIJr1wedvrr9MADD9izLv/880/Gjx/PG2+8weWXX86SJUs8JqAvo1AoXHLXyFYLp//4mh2fP8/WrVu9pp5Sc1ARo2Ku8oBjNFuZ/u1JDmYZMFlkFu8p5I1VmRSWWXhtRQY3fXmC+RtysFrrd9HN5QagL8fD1Bd7wWh3x6Rpz5uvleLt+tNoNDjEyjYV+Ye3sXfeo+x8707OrP0OWXbuoZAVNh+BPcGqFfDGykx+22trsH4wy8D0b09iNFv5ZVc+t3x1gicXnyW9oOUnMjUUs9nMRx99xP33308RRbYkAVcNtAoUIPeVsfa2snHTRqbeNJUdO3Y0pbhNglNP2oEDB+wB/e+99x5ffvkliYmJ5ObmMnXqVC699FKPCemrKJVKZGvtTzPG4jwUSjUqXSDZ21eStfV3AMzYlnuvueYaIiIiPCCtd2E30qrYW8dzjZzTO17PradKuPP7U6w5bKtJt+pgEfmlZu67qO7SAiar41wtGbvR5HqnGdfQnDdfK6W16E+LsQxzaRHakChMxfkcW/SevcRQ5qbFaEOiiOx5UTNL6V40Gg1Fcv2WTLaecuzJfE5vYd6GHF5enmnf9ueRYlbemdrgWmkma8t8+Dx+/DgvvPgCBw8cxNrWitxbbmBUfRUkkFNkzOFmcjblcM899zBlyhRuvvlmr6ml5tSTZrVa7cVYJUmyB3aGh4e3aDdqU+Ln5wey1fYPMJcWYy1XYFaziSML3mH3nHvZ+f7dnF7zLUWn9zucX1ZWxubNmz0utzdQEXRusFQqqrbhGiIDHJ/+u8Xp7AZaBQt21s9AqBi7NQS45+fnI2mkBkahuoDwpAHerz81Gg3U8QB5PlaTEYuh0sg4t3M1u96/mz0fPsD+L56h4Ng/DjUgAYpOOy+xIZXP7y0/hg3Fz8+PMmv9vJJ9kxzjAKMCVfx1zFF/Hcsx8s+ZhpdfqtBrvn5dKzCbzXz++efcdPNNHDpxCMsgC3L/JjDQqhIO5ovNWNpZ+Prrr7np5pvYtWtXE07QcJx+zBkzZvCvf/2L66+/nj59+nDPPfcwatQoNmzYwLBhwzwpo89SkWpvKSng2O/zKDy2E6XWn8QLr8NqNlFweJvtQKuFrM2/EdlzpMP5arWaXr16eVhq76Di2pWZK400jUrBnGvb8MTiipi0YB4cFc3iPQUO/Ttjg+t395aZJbQatVcUF3U3RUVFNiPN3agr52vNeLv+DAwMQsqrf2B55ubfSP9rAVazkbCOF5Bw4RROrfwC2WJbmivJOEZRWBxICvtDKUBAbDvng1ps2esV3Rh8lcDAQErM9bu3HhwVQ47ebI9Je2F8PD/+kw9UxuVJEkQHNdwCqZAlMDCwjiO9n7179/LKq69w/NhxrInl3jN3PVOrbMuflkQLZ7aeYcaMGUyaNInbbrutWa+l02/CuHHj6Nq1K9999x3Hjx/HYrGwfft2xo8f7xVKxheo+MNmbvqNwmM7AVvfu5PLPyG865Bqx/tHJRHeZTB5+zagVquYPXt2qw2+rrh2+vOUX//kAJbekeqw7fExsTyxOB2TRSZEp+Thi+t3zfRmqUUosvpQVFSErPZAIooECq2CwsJC98/lxXi7/oyKjOBo9vF6HVt67gxn1nxrf593YCPqwDC7gVaBqTiP5LH/4czab7GU6QnvOoSoXs7LKymMNq+cr4dzhISEUGiwJSLVlZQe5KfkvWvaOGyLDVbz55FiTuebkCS4Y2gUiaENj0soNNqECA0NbfAYzU1RUREffvghCxcuBB1YBlsgwUOTx4B5jBlpt8TCXxay9s+13HP3PVx00UXN0sGhVnM9OTm5xp6agvpR0dKpNPuUw3bZakEb5hjtKCnVBHfoRWSvkaR064M2/R+mTp3qKVG9joosm2JT3TfFlL7hjOwYxOFzBnol+OOvqd+aXrFJajVtt4qKirCqPZRmrq7sW9ua8Wb9GRsbi/Kf3fU6tvTcqWrbzGXFqIPCMRVVtvQL6dCLiG5DbQ+gshWpjsQEyVCEJEn2FoG+SkREBEYLlFpsDc5dJSlMw+q7O7L9dAmxQWrahDcucDTfoECjVvvkA6gsy6xatYp33n2H/Px8rB2syN1ku4feY6hA7iVjaWMhb1sezzzzDIsXL2bmzJker+nn9Nfs5ZdfZuvWrZ6UpcVRETgcEt/eYbvSL4CIHiNoc8l/8I/rQFCbLqRcORNNkO14hcVAcEhoq1iGc4afnx9+Wg2FxvoZXNFBaga3C6y3gQZQaFISHuFaWxxfpai4yGOKTlbLPl9WobF4u/5MTExENpaCqe7MyqDETkhKxy9PcNtudJh8HyEpfdBFtyF+2FVE9x0D2GLw6jLQABRl+URGRft87FRMjC1JKaes4QGfaqXEBckBjTbQAM6VKYiOivS5vp1nzpxh5syZPPvss+RL+VhGWWzLm82Z/xAOlpEWrL2sbNmxhRtvvJHPP//co60Endr9CxcuZPPmzeTl5XHppZcyYcIEunTp4jHBWgIVbvyY1O4YUZG3/2/UgWEoNDp2z74Lpdaf+GFXEdnjQofzJKOeiBbaKqu+SJJEWFgY+YbqHpkj5wzM35BDmdnKDf3C6Z3YsJiWAqOS5FZynfV6PbK/Z+ruWZXWVm+kebv+7NChAwCKkhysIbWvI6kDQ+kw+R7Orv8Zi6EE/5i2nF79NeaSIsLSLiDtuidr7OlZF6rSXDr27tog+b2JuDjbqkhmiZKkwLq91WUmK59szGHX2VIGtwvkur5hKBqYyVkTWWUq4jvVVXrfezCbzXzzzTfMnz8fM2asvcsL03qLjakAOVXGnGjGusPKhx9+yLLly3j4oYfp1q2b26d3aqTFxsby008/cfz4cRYvXsyDDz6IxWJhwoQJjB8/nnbtagkIFQAQGWl7mlGYSogfegXxQ68ga9syTq+yNVg2lxZxcvknBLXpgja00uWvMpUQF5fmbNhWQ3RMDLknzjhsy9WbueKjo+SX2jLDFu4s4JdpHegc61o0qVWG3DKIinJzx3EvobS0FDxVu1QN+pLWbaR5u/5MS7PpF2Vxdp1GGkBw2+4Et+2OSV/A7rn327M48/ZvwC8inrhBk1wTwGyEkny7HL5MUpLNIEovqZ8n7aGFZ1i4y5b9/OueQs4UmHjo4rpLBtUHqwwZJUr6JvmGkXbkyBFeePEFjhw+gpxgq1lG87e2rRkdWAdZ4Syc3HGSGTNmcNVVV3Hrrbe6tUKA029Vhau0bdu2zJgxg8WLF/Pf//4Xg8Fgb+MkqB2VSkVYeDhSFW9QScYxx4NkmZLMEw6bJGOxz8dpNAUxMbHkGByfI5YfKLIbaABGi8zCXfkuj11klDBbaTXX2VBmAA+tnssqmZKSkroPbMF4u/4MCQkhMSkJRXFm3QdXoTT7VLUyG9V0Wj1Qls/bvXt3l8/1NgIDA4mKCOd0cd03WJnJyq97HMvT/LAjr8lkOVemoMwsN/tDQF1YLBa++OILbrnlFo6eOYplkAXrYC820KoSD+bRZiztLXz//fdMvWkqe/fuddt0To20mloSderUiZkzZ7J8+XK3CdTSSEhIQGmoLEcQmNTZYb+kVBGQkFK5wWxANhl8uuFwUxEbG0tOmWPXgfPrpAFolBL3/niKS2Yf5sWlGZSZ6l5yyCq1ffVbw3W2WCyYzeamrStUG8qW3zewLnxBf/bp3Rt1cQa40LfQP7YdCrVjDFlgUieX51YWnEWpUtG1q+8vdwKkdOzICX3dS75qpUTweW3sIgJUzFmXzbg5h7nlqxPsz2x4B4YTRbaxK5azvZFz585x7333MnfuXExxJsyjzZDY3FK5iBrkPjKWERbS89O5Y8YdfPnll27pAerUSPvyyy+bfLLWSFJiIkpjZTmCiG7DiB14GaqAEPwiEmg3cQaawDD7fkWZ7SkrIcFT+cbeS0JCArIM2aWVX9MLU4MYkVKZtdQ5xo8/DhXx884C9mWWMfevc7yyvG7vQGap0j5HS8duMLngSTu37Ry7/ruLPf/bQ/6BfNcmVILRaHTtnBaGL+jP/v37I5tNLnnTVH4BtL/8HnTRbVD5BxPd9xKi+4x2eW510Rl6dO/RYgpJp6V14kyxRJm59uOUConHxsSiLFdpfmqJvkn+vLw8kz3pZSw/UMSNnx/HYG7Yj/3RQiVKpYKUlJS6D24Gdu7cyU0338TO3Tux9rciD5Td0qquNLuUA/MP8M9r/3Bq6SnkerYKdJloWxFcc5yZDz74gMcef6zJVxGcPlsHBAQgyzI7d+4kMzPTnirdo0cPn8saaU4SExORDSW2GAyVxlZ9fOiVxA+9ssbjFaU2Iy3JR2IK3EnFNcgoURIXYFNaSoXEZze2ZdupEsrMVlIitfR/w7Gq+YqDhTwzrvaGbul6BQqFwh7025KxZyLVM/ms8EghBz856PC+9xO98Yuo5w+qQhhpvqA/+/Xrh1KpRJl/Emtw/e+D4OSuBP/r+QbPKxmKQJ/LoEFTGjyGt9G1a1dkGY4UqugaXruldk3vMIa2D2BvRhl9Ev2596fTDvuziszsPFtK/zauN14/VKAmNSXFKzNmV6xYwYsvvojV32rruemmGFnZKrNvzj7Ksm0eSf1pPZJCInG0m9x1GpAHylgPW/nrr7+4ffrtvPnGm00W7+xUba9bt44xY8Ywa9Ys1qxZw+rVq3n33XcZM2YM69ata5LJWwPJycmALd28Pkhl+SgUilbh4amLimt3Rl/9a9onyZ/B7QIJ91cRFej4rNEhsm4FdUavJCE+rsX1t6sJV4203N25Du9li0z+3vz6T6jAK1ofNSe+oD8DAgLo3bsPmvyTtkqsHkKZZ4vBHTKkekFvX6Vr164oJIn9efWLKYgP0XBxWjDhAapq+kqlgOQw17NljRabkdijZy+Xz3U3ixYt4rnnnsMSbnGrgQZQklFiN9AqyNvVdHF/NSLZMkAtQy2cOH2C6XdMJz09vUmGdvqNevHFF5k/fz6JiY7W56lTp5g2bRpLlixpEgFaOnYjrTQfa2DdQeqK0nzi4uNbhfFQF0FBQUSEh3Gq2Hl8k0op8cpl8cz8+Qz5pRbaRWh4YmzdHQdOl2jo2Nl74zaaErvBVE8Hji66evSuLsaFiF4FWMyt20jzFf05YsRwtmzZjFSai+zvmcr/6rzjtG3XrkWtFgQGBpKamsKerH1ciWsxZXcOj2L76RK2ny5Fp5Z4dHQs0UGu6//DBSpMFujdu7fL57qTlStX8vrrryPHybbsSDcnMGnDtCg0CqzGyiVjvxgPLavHgnm4mew/s7nn3nv4YM4HhIWF1X1eLTh9trZYLDW2JIqJibEFIQvqRXx8PCqVCqm0fpa8ylBAey/PzPEkKakdOamvXWGNTA3i6UvjuKFfOK9MjCc1qvYbsswCmXq8Nm6jqbEHsdfTSIvqH0V4z/L6cQqIHRZLSEcXHn2lmgPnWxO+oj+HDx+OJEmoclzP0GwIkkGPojCDkRdd5JH5PMkFAwZyuECFvh5dUqoSEaDi5YkJTB0QziMXx3Btn4b9qO/MUaNSKr3KSDtw4AAvvvQiROERAw1ApVPR7qp2KLQ288Y/3p8249rUcVYTEg7moWYyszJ54sknGn2/O/WkXXnllVx11VWMGzfOHreTnp7Ob7/9xlVXXdWoSVsTKpWKpKQ2HMrLr/tgqxlKC2jfvn3dx7YSUlNT2bJpI0YLaJzc4E8sPsuXW2xG8JdbcnltUkKtiu5kkRIZ786AakrsGUf1/O1QqBR0+k8nDHkGJJWEJqhhVdCtVisKRcOrsPsyvqI/w8LC6NO3L1t3H8SU2Lfu5pONRJl7FIBRo5z39PRVBg4cyOeff86uHBUDY+tfkf7vY8X832fH7Vnsv+8v4puprj+ob8/R0r1Hd69pWG8wGHjm2WewqC1YBlk8VgIIIGZgDJG9IzEWGtFFNUNdjwiw9LWwa9MuvvnmG/7v//6vwUM51aC33XYbb775JgA7duxg+/btALzxxhvcdtttDZ6wNdKhQ3vU9YhJU5QWgCzTtm1bt8vkK6SlpWGR4ZSTGkRFZRa+3ebopfz473O1jnms0PZs0rlz51qPayk0NFBdG6ZtsIHWmHlbAr6kP8eMHg1lhSiKs90+lybnCKmpHVvUUmcFXbp0ITQkmC3Zrt0zn27KdSgz9PcxPXvSS10aI6NEwZliiWHDhrt0njv58ccfOXP6DOa+ZrdkcNaFUqtsHgOtHDlZRk6QmT9/PtnZDb+3ao1y7NChQ6vxNriT9u3bs2LFCnuGpzOkUlvAtrjmlXTqZKvBdKRARYeQ6nFOCsmW8WmukmKtVtZuHBwpVBIRFkpkZOvo22n3ZnlqBVIu79/Yio008B39OWzYMN54401MOYcxBrmvuLNUmoekP8fYsS0nq7MqSqWSYcNHsGzJrxgserT19Bypa2gJVZcOO5+NmbbflWHDhrl0nrswGo18/c3XyLEy1B0i3GKx9rBi+t3EDz/8wPTp0xs0hlNP2tq1a+2vi4qKePzxx5k4cSIzZ87k3LnaPRUCRyo8Y4o64tIUJXkolcpqwcatmejoaCLCQjlcWLPGC9AquXlgZcCzQoI7htWe+nykSEuXbr5f6by+KJXl167p6yzWjAxSE/Yi9EV8SX8GBgYydOgQNLnHXCps6yqqc4eRJKlFLnVWMGrUKAxmme3Z9Q/8v2VwBDp15f0ytlMQHaPrH+guy/B3ph/dunaxN3tvbjZv3kxBfgHWFE8pHS8lEOQ4md+X/t7gOF2nRtrbb79tf/3KK68QGRnJnDlz6N69O0899VSDJmutVMSY1WmkleaRmJSESuWp0vDejyRJdO3eg0OFzv3lj4yO5at/t+WpS2L5fXoK47s6D3IvMEpk6vFIY1xvwf598pQnzQpqVevOTvY1/Tl27FhkUynKgtN1H9wQZBlN7lH69u1HRIRnskibg549exIZEc6f6fVf3+uZ4M/KO1N55tI4PryuDe9f61qQ+4kiJaeLJcZecqmr4rqNbdu2Iakk8A6bsVmR42XycvM4ceJE3QfXQL2ienfv3s19991HQkICU6dO5cyZM3WfJLATGxuLWqNBUZpf63FqkdlZIz169CC7BHLLnHtnhrQP5D+DIkmrI9X6YL7NYGkJPQPri72ci6eqYlhApRYPGhX4gv684IILCAwKRnXusFvGVxRlQlkRl1wy1i3jewtKpZKxl1zKrlx1rfrqfBJCNdw0MIIxnYJRuuiFXnNWg1qt4iIvypg9efIkcpBc79qMLRk51PZ0fPLkyQad7/QS5uTkMH/+fObNm0dxcbGDq84d/alaMgqFguTkZKTajDSLGbm00Osb4zYHFQbVgfzG//AfyFOhUavp2LFjo8fyFezVxz1opHljxXNP4mv6U6VScfGokagLToKl6btFqHIOo9FqGTp0aJOP7W2MHz8eq4xL3rSGYrDAX5k6Roy4kODgYLfPV1/0ej2yunWX4bFT/oys1+sbdLpTI+2aa65Br9dTUlLC5MmTycuzLdVlZ2e3mqy4pqRd27aoDflO91f07KwofiuoJDU1FZ2flv15jV9C21+goWu3rmg0Dc9a9DVUKpVtydODRlpL6cnYUHxRf44ePRrZYkaV27BlGadYrWjyjjNs6FCvKQ/hThITE+nbtw+rzuqwuNke35ChQW+Sueyyy9w7kYv4+/sjuVgvrsVSXo1Fp2tYpqlT18Sdd95Z4/aoqChee+21Bk3WmmnTpg1y2XKwmEBZ3diQykt0tGnjwaJ7PoJKZWt1snf3BooNRfywI5/8EgsTu4fYW6qYLTI5ejMxwc4NOb1J4kShgot6eU+xR0+h9dNiNHmmn6ZkkvAPavk/xrXhi/qzW7duREZFk5l7FHNUapONqyw8g2wq4+KLL26yMb2dyZOv4Imt29h2Tk3/aMeaaX8cKmLTcT29k/wZ06nS+5VdbCZIq8BPXb81QlmGZad1tGubTM+ePZtU/saSlJTE5m2bbclKrXzJUyqwGasNLTtT6/rRkSNHWLlyJVlZWYAt027UqFE+kVbubVQYX4qyQqwB1QNnFaUFSJIkMjud0KdPHzZs2MDkj45xMMvWduX99dn89J/25OjN3P/zGbKLzaRFa5k7pQ1tI6ovNezLUyGXj9XaCAgIoMhU5JnJzBAUGOSZubwYX9OfkiRx8aiRfPPdd2A2gKppluuUOUfR+QfQv3//JhnPFxg8eDAx0VEsPWV2MNLe/zObV1Zk2t/fOTyKWwdFcPt3p/j7mJ4grYInxsYypW94nXMcyFdxokjBg7df43Xlbnr16sWPP/4I2YjkgXQIDglucCiTUxt37ty53H///YAtJqgiLuj+++9n7ty5DZqsNVNhfEllhTXuV5QVEhEZ2epjeZzRp08f8vPz7QYaQJlJ5ovNucxcYDPQAA5kGXju94wax9ibp8JPq6FLly4ekdmbCAwM9Njyg8KsIDAw0CNzeSu+qj8vvPBCsFpR5TXRkqfViqbgJMOGDml1IQZXXnU1+/NUHKtSPujjDTkOx83fkMO7a7L4+5gtXqnIYOXJxemcK667ldCSk34EBwUyevTophW+CRg4cCABgQEoDrdyN1oJKM8qGX3x6AZ3X3HqSfvxxx/59ddfqzX6njp1KhMmTGDatGkNmrC1kpCQAIDCUFhjaJDCWERSivCiOSMlJYUA/+pr+hYZsoocFdqBrJobHO/J09KjR89W2bw+LDQMKd81I81cakZSSCjrW5WzHMkoeVUQc3Pgq/qzc+fORERGkpl3AnNU45NrFEXpyCYDI0aMaALpfIsJEybwyfx5LD5h5M7uNiPs/CK1KqXEwSyDwzajReZYjoHIQOcLXel6Bduy1fzfjZO9Mv5Tq9Vy7TXXMm/ePJs3rfbSlW7BarFiKbWgDmw+fS/tllBICq655poGj+HUtJMkye6mr0p2drbXuVZ9AX9/f4KCgpEMNS85KY1FxMfHe1gq30GpVHLx6DFEhFX++Af7Kbh9SCTd4hyV1IWp1Zfa8gwSZ4ol+vbr53ZZvZGQkBAUpnrGulhlDn91mE2PbmLzY5s5+asLqeMyWMushIS40JC9BeKr+lOSJIYPG4a68Iytl3AjUeWdQK3RtKqlzgoCAwO5bNLlbMrSkFVqu/fuGu5ordw5LIqLOjrqq4gAJd3jaw8y/+2kHyq1iiuvvLJphW5Crr32WqJjolFtVdmD5z3Fue3n2PLkFjY/tpld/92Fscgz8bgOnAHFCQXXX3+9vX9vQ3Bqqj/22GNMnTqV5ORk+wRnz57l5MmTPPnkkw2esDUTFxdHXmYNRprVjGwoITa2FffPqAf9+/dndZfujIs4hWwxcUmXEKICVXwwpQ0vLM1gX0YZI1ICeWR09eu4O8f2NNW3b19Pi+0VhIaGQs0Oxmqc23qOrA02A8NqtXJ62WlCO4cS3KEe3jETYC2frxXjy/pz0KBB/PzzzygL07GENqLHpiyjLjhN3z59vdLb4wmuuuoqfvj+e5ac0PLvTqVc3y+c7vE6Np3Q0yvRn75J/litMsUGKwt3FRAfouLhi2NrTR7IN0isS9dy6YRxhIfXHbvWXOh0Op568inuvuduFBsUWIdYPZJEYDFYOPLVESwG25pV0dEiTi0+RYcpHowFzQfVZhUpHVP497//3aihnBppw4cPZ+nSpezcuZPMzExkWSY2Npbu3btXtpkRuERMTDQHTmVydt0PnNu1FpVfAPHDriYs0fbliY52X9+8lkDfvn1RKBTEREcxOqlyiSAxVMOcOqp078lVERocREpKirvF9ErCw8OxGq22Mhx13L76s9Xr+ejP6OtnpBkq52vN+LL+7N27N2q1BlPB6RqNtMLjuzmz5ltMxXmEdxlMwohrkRTVP5NkKISyQgYOHOAJsb2S6OhoRo8Zw8plvzO5fRnBGpnu8ToHT5lCIXHPhdHcc2H99P/SU1rMssSUKd7fA7VHjx7MvH8mr7/+OtImCfkC9xe4LcspsxtoFejPNKxGWYMoBNWfKsKCw3j5pZcbHYtZ6+VSKBQkJiaSmJhImzZtSExM9HoF481ER0eTdfIwGRsWYdYXUJZzlmOL3sOUl27fL3BOQkICsTHR7Ml1raitLMPufC19+vVvcPCmr2M3murhTQvtFOq4QQGhaaE1HVqdUtt/Lbn1T33xVf2p1Wrp0bMH6sKz1faZy/QcXfgupdknMZcWkbV1KVnbltc4jrLAdn5rXOqsynXXXYfRIrPsVOOTwkrMsOKMPyNGjPCZSgATJ05k+vTpKE4pUPylgMavoteKf4w/mjBHwyi0c6h7J60gB1SrVYToQnjnv+80yW+601+7ffv28fTTT1NUVERsbCyyLJORkUFwcDBPP/00Xbt2bfTkrY2IiAjycx2ze2SLGf2ZQ4QgftjqQpIk+vW/gFVLs7HKeurbPeWsXkF+GfRrpfFoYKvPBdiMqIDajw1NC6XDlA6cXXMWhUpB4phEdDH1K8QolUqO87VSfF1/9u3Th61btoCpFNSVf/uSjGNYTY6B7sWn9hPT75JqYygL0wmPiPAZY8JdtG3blqFDh7J803omJJfh14jGKStPayk1ydxwww1NJ6AHuO666/D39+ett95C8YcC82BznXqooUhKic63debEghOUZpcS0TOCxDHu/w5KxyWU25TERMfw1ptvNdn33unX5ZFHHuG5556rViRvx44dPProo/zyyy9NIkBrIiIigsDAQMeAYkkiIDgU8jKFkVYP+vbty6+//sqxQiUdQupXQn93buuOR4NKo0kqlZDr0Wk9ZnAMMYMbUOCo3JMWGRnp+rktCF/XnxVyK4sysYS3tW/XRbVBUqqQLZXuEP/YGuo/yTJqfSZ9hg3y6kQJT3HDDTewbt06Vp/VckkbQ90n1IDJCktP+dO3bx/S0tKaWEL3M2nSJGJjY3n6macpXVGKuZ8ZEtwzV0B8AF3u8FCpJTNIOyQUxxT07N2T5559rkljcp2u/ZSWltZYxbhXr16UlpY2mQCtidDQUOLi4ojsdAGSQolKF0SbMTeh8/NDqVS1+tpS9aGiEO2e3PqnVe/JVREfG9OoDBtfJyam3OAqcfNEJRAQGNAq2v/Uhq/rz44dO6JUqVAUZzpsVwcEk3zpragCQkBSENZpIDH9Lq12vmQsRjbo6datm6dE9mq6du1Kj+7d+f2UP+YGtor6K11DvgGuv963vGhVGTBgAPM+nkdK2xSUfymRtkgez/xsUnJBtVKF4piCG264gbfefKvJk6ZqTRyYNm0al19+uT3rMCMjgwULFjBs2LAmFaK1EBoaikKhoN2QiSRdehsoFEiSAunonwQGBYknznoQFhZG+3Zt2Zt3iMvqUcDZKsP+Ai0jx7buuJiAgAACAgMo0ru364Ckl1q1MVyBr+tPrVZLSkoKe8+eq/YbGt5pIGFpA5CtFhTKmn9CFMXnALy2T2lzcN311/Poo4+yOUvNoFjXLBNZht9O+dOhfXufD9uIj49nzvtz+Oijj/jmm28gG8x9zL7VmcAC0l4JxQEFYeFhPPn2k25bqXFqpD3xxBOsWbPG3tZElmViYmK44YYbWmVhwqagosCnZDYgVVFuktlASEjrLv7pCr379GXRghOYraCqIw/gRJGSEpPcKltBnU9sbCzFxcX1Wu5sKIoSBfEdRb2/lqA/O3fqxMHDv9kshPMeICVJctBh56MoOYdCqaR9+/buFtNnGDRoEEkJCfx+6hSDYgtcOndXroozxRKP3T2lRTzMq9Vqpk+fztChQ3nxpRc5u/Ys1rZW5J4yeHtjimxQbVMhF8pccskl3HXXXQQFua8NntO7zGAw0L1792oKJScnB4PBINoXNQD7H9LiGJMgWQwEB4l4tPrSu3dvfvzxR44WKukYWntc2r4821e8V69eHpDMu2mT1IZjW49hpYHrLXUhA8WV3TVaMy1Bf3bo0AHZbLQtXWpd+xFSlOSSlJTkE5/TUygUCq68+mr++9//crhASUo9Y2oBlp7yIyw0hJEjR7pRQs/TvXt3Pv3kU+bPn88333yDnCFj6WFBbiODt9miBpB22WLPomKieOiph7jgggvcPq1TP8QLL7zAli1bqm1fv349L730kluFaqlUxOlIFkdXt8JqFvFoLtCjRw8A9ufXnSZ1IF9FQnxcqw9kB9syg7XYirtsNEpt3QqEkdYy9GdFQ2hFSZ7L56rLCkjx0kbyzckll1yCv86PFS6U48gqUbDznJpJl09ukf1PtVott99+Ox999BFp7dJQbFKgXKsE90Zm1B/ZlrmpWqZCdULFlClT+Pyzzz1ioEEtRtrWrVsZM2ZMte2XXXZZjcpHUDcqlQq1RgMWxxYVCoux1Qdau0JoaCjJbZI4kFd78oAsw8ECLT179faQZN5NUlKSzUBzV/JAUZV5WjktQX8mJycDIJW5tjSH1YxcVmg/X1CJv78/Y8ZewsYsLcWm+rmK/jirQZIkJkyY4GbpmpeUlBTmvD+HmTNn4l/sj2q5CmmPRI3Nrj1FISjXKlFsVtCpfSc+/vhj7rjjDnS6+pUkagqcGmmy7DxuxWp116N4y8dP64dkOa+an9Xs0T96S6B7j54cLtJgrSW8Kr1EQZFRFhlm5diNJzc9oUpFkuM8rZiWoD9DQkIICAhEUVbo0nlSme0LJjyqNTNx4kRMVvg7o26vmFWGdRk6Bgwc2CqKnSsUCiZNmsSXX37JqItGodirQLVcBdXb4LoXC0h7JFQrVPjr/XnggQd4f/b7dGgG77BTIy0iIoKdO3dW275z585W3/KlMWj9/MB6XmaPxdxqe9s1lG7duqE3yqTrnWcOHCpQ2Y8VVPGMFLop2KMQdP46Ue+PlqM/4+LjkAyuWfWK8uPj40UCSU2kpqaS0qE9f6bXrfN356rIK4Nx48Z5QDLvISIigqeeeoo333yTmMAYlGs8WK4jB1QrVCj2Khh54Ui++vIrLrvssmbrVuM0qOehhx7i3nvvZfLkyfbq2Lt372bBggW8/fbbHhOwpaHVapFKHP23stUsAmxdpEsXW6HCI4UqEgKNNR5ztECFv05Hmza19/VsLQQHBxMaFkpuYa5bxpeKJNq2bdsiss8aS0vRn3GxsRw+sxtXyq9KxmKgSm0+QTXGjL2E2bNnk1GiINbfuWf1rwwNgQH+DBo0yIPSeQ/9+/fns08/sycWkImtCK47vlrl3jPFQQURkRE89NpDDBw40A0TuYZT07BHjx589913yLLMzz//zIIFCwB49dVX7a8bSn5+PjfddBNjxozhpptuoqCgesxDeno6N954I5deeinjx4/n008/bdSc3oJWqwFrFSNNlsFqaZEBoe4kKSkJf50fxwqd90I8WqSmU+dOrbZfZ020a9cORaF7roeyUEm7tvUoXtcKcKf+9CTR0dFIRteaU0tGPQql0qc8hp7moosuAmBzlvO4WrMVtp3zY9jwEa3698HPz4/p06fz/vvvkxCRgHKtEmlHE8eqFYJqlQrFAQXjx43n888+9woDDeposB4ZGcndd9/N9OnTSUxMZMGCBbz77ruNXpedO3cugwYNYtmyZQwaNIi5c+dWO0apVPLII4+wZMkSvv32W7766isOHz7cqHm9AY1GA3KVJ6fy12p1/SvoC2zfj9TUjhwvrvm6ma1wqlhBWlonD0vm3bRv19623NnUpdLKwFpmFXWxquAu/elJIiIikM1GOD+OthYkUykhIaHi4agWYmJiSOuYytZs5yso+/NUlJhknyh+7Am6dOnC/HnzmTx5MopDCpSrleDa80ONSCckVCtVBFoDeeWVV3j44YcJCHBTY9EG4HS589ixYyxevJjFixcTGhrKuHHjkGWZzz//vNGTrly50j7O5Zdfzo033siDDz7ocEx0dLQ9UDIwMJD27duTmZlJSkpKo+dvTjRqNZJc1ZNmM9JUqkZ03W2lpKSm8uueXVhlqjVbP6tXYrba4j8EldhqX8k25daUVV8KKscXuFd/epKwsDAAJHMpsrJ+tdIkUynhkWHuFKtFMGjwED779BBFRokgTfWnpn9y1KjVKp/vMNCUaLVa7rvvPvr168cLL75A2coyzIPMENWAwWSQdtqWN7t178azzz7rlaWanD7qXHrppWzYsIE5c+bw9ddfc+ONNzbZk1FOTo7dAIuOjiY3t/YYmdOnT7Nv374ae+H5Gkql0rbEWUG5kaZUOl+2E9RMhw4dMFhkskqrfy9PFivtxwgqsXu6XKyqUBdSgeQ4fivHnfrTk1T0IZRMZfU+R2EuIzxMGGl10b9/f6xyZcHt89mdp6FH9x4iqawGhg0bxkcffkR8VDyqP1Vw2sUBLKDYoEBxUMHkyZN55513vNJAg1o8abNmzWLx4sX861//YtiwYYwfP77WtPLzmTp1KufOnau2/d5773VJQL1ez913381jjz3WIgq+KpVKpBqWO4WR5joVxTbPFCurBd+e0StQKhUkJiY2h2heS/v27W0tffIl5IQmXPPMh9DwULvnpbXTWP3pLVRtZVdflBYjISEh7hKpxdC5c2f8tBr25am4IMYxbbHYJHGqSMGY3qLGozOSkpL4YM4HPPzIw+zZsAfrBVZbp4K6sIDybyWkw4wZM7j22mvdL2wjcGqkjR49mtGjR1NSUsKKFSv45JNPyMnJ4emnn2b06NEMHTq01oE/+eQTp/siIiLIysoiOjqarKwspwGmJpOJu+++m4kTJ9ZYGNIXUSqVSFUCgqRyxS2WO12nbdu2AJwpUXB+a9t0vZLEhHgR63cefn5+JCQmcCr/VJP28FQWKOnYsWOTjefrNFZ/egvOWtnVisXYIh6o3Y1KpaJT584cObYNKHXYd6TA9tAuygfVTnBwMG+9+RYPPfwQ/2z6B4vaAnG1nCCDtEmCdHjggQe47LLLPCZrQ6nT/+7v789ll13GBx98wJo1a+jcuXONgf6uMHLkSHuG04IFCxg1alS1Y2RZ5vHHH6d9+/bcdNNNjZrPm7AteVT9cZSrbBe4QkBAAGGhIWSWVPdCZpSpSWrT1vNC+QCd0jqhLGhCz60F5EKZtLS0phuzheAO/elJKgKoJXPNZW6qIcvIJoMw0upJ585dOFlki5+tyvEi20O7uKfqRqfT8eorr9KhQwdUG1VQS+1laa+E4rSC6dOn+4SBBvUw0qoSGhrKlClT+Oyzzxo16bRp01i/fj1jxoxh/fr1TJs2DYDMzExuvfVWwNZWZeHChWzYsIFJkyYxadIk1qxZ06h5vQGFQuHYN1YWRlpjSExKIuM8I80qQ2aJJCqeOyE1NRWr3opLxa9qowAQSRp10lT605PYs9zOL8DtDNkCslW0uasnHTp0wGyFzBJH/X+qWElsTJQwduuJv78/r77yKoG6QJTrbUuZZDr+kw5KKPYqGDt2LFOmTGlegV2gWdbYwsLCaqx7FhMTw4cffghAv379OHDggKdFczuSJDkmDpR70kQB0IYRH5/AlsN7HLYVGCVMFlHx3Bn2Zck8ILbx40n5kuO4ghaDTqezxTBa6mmklR8njLT6UdEF5GyJkoTASndaRqmK5M4iCccVoqOjefSRR3nsscdQrqt5pSA2Lpb77rvPp35vRSCUh7F9Oapmd1bdLnCV2NhYcstkLFZQlj+M5pQp7PsE1akwpqQ8CTm2CeLS8sA/wJ+4uNqCQQS+iCRJaLRajPU00iRhpLlExYPk+RnqWaVK+oiVAJcZOnQoX375JXl5eTXub9++vc99N4WR5mFkWebIvl1krPwdpV8AiYMm4o8w0hpKVFQUsgz5RokIP5vBkVtupLWGhsQNISgoiNi4WNLz0pskeUCZryQtLU18h1sofn5+FFcx0vIPb+P0H19h0ucT3nkwSRf/C4Wy/KekfFlUlI2oH0FBQej8tOSWVZY4KTVDiUkW+quBJCUlkZSU1NxiNBkiEMrD7Nixg7MnjmI1GTAV5XJs+eeUlZWJmLQGEhVlq2KYZ6i8frnlr7217o030KVzF5T5TZA8YAHyoXOnzo0fS+CV+Pn5IVltHQfMJUUc+3U2xoJsZLOJnF1ryNqy1H6sVN6ZQKfTNYusvkhEeDgFxkr9VfE6IiKiuUQSeBHCMvAwZ86ccdwgWykqKmoeYVoAFeVbCgxVlZyEUqmw13gSVCctLa1pkgcKQLbKdOok2m+1VHR+Oig30kqyTiCbHZc+9WcPVb4pP0540upPcGgoxaZKL3TFa6G/BCCMNI9TzQ2rUBAUFCQ8aQ2koiJ6gbFSyRUaFYQGB4vlt1qwG1W1N/uoEylXchxP0OLQ+evsHjL/mLZIKsfagwEJlVm9kjDSXCYwMIhSS6X+LzVL5dtFZqdAGGkep1+/fiS0aYdC7YcmOJJ2F/+fbTlBGBQNosJIKzY7PomGlG8X1ExFDJmU18jvXR4EhwQTExPTNIIJvA6dnx+SbDO+VLpA2k24A21oDAqVhsgeFxLdd2zlwRZhpLmKVqvFaK38KS6z2O5JcQ0FIBIHPI5SqaRDp67EXPMMAFJJHuw6JjxpDUSr1aJRqyk2VV6/YpOCoGDRlqY2/P39SUxK5GTuyUYlDyjzlHTp3EU8ZLRgdDodinIPGUBoSh9CU/rUeKwkEgdcRqPRYLJW3j+W8kocoluKAIQnzeM46zggfuQaTmCAv32JAKDEoqxsZyNwij15oKE2mtnWaUAsdbZsdDqdfRmzTiqWRX2szEFzolAoqNpwoOK16OcsAGGkeZzzi9lKouNAowkICKCkipFWapEqK6ULnNK5c2espdbz2wbWnzxAto0jaLnodLp6dxyo8KSJ7M76Y/tNqL5dlpuut67AdxGWgYdRKpUgV31uEkZaY9EFBFBWxUgrM0viSb4eNDZ5QCQNtA78/f2hvr07LSZUajUqlYikqS9WqxVFlYWUipfCSBOAMNI8js1Iq3LzWW0Gm1BqDUen88dgqXxvMMsiJqYedOjQAYVSYTe2XCYPIqMiCQsLa1rBBF5FQEAAssVs11W1IVmM4gHJRSwWi4ORpix/bTbXc4lZ0KIRRpqHUalU53nSbK9F/EHD8fPzw2i1XT+rDCarCFyuD1qtlvbt2zc4w1OVp6Jrl65NLJXA27CHDljq9qbZjDQRauAKJpMJlVT54F7R3k4YaQIQRprHsRlpVdw+5U+nIpOn4Wg0GoyyzdAwltu/Wq22GSXyHbp07oIiX+F68oABrMVWsdTZCqio1yXVw0jDbCRYJO24hMlkQq2ofHBXK2w3o9FYzyVmQYtGGGkeRqVSIVdZNpDKDTbhSWs4Wq3WnsJuKq8xJIze+tGpUydkowzFLp6YV3m+oGVTkSktmetuT6GwGAkJEZXyXcFgMKCSqhpptv+FkSYAYaR5HLVaDVZLZVxa+dKnRqNpRql8m6p1hszll1V40upHhZHl6pJnxfEdO3ZscpkE3kVIiK3mYH2MNKXFYD9eUD+MhjI0VX6JhSdNUBVhpHkYuzFWEZdmtXnShOen4ajVaszll9NsFZ40V2jbtq3tWrmY4SnlScTFx4l6dK0Au9FlLqv7YHOZMNJcxGAoQ6OsjDeoMNgMhsY21hW0BISR5mHsRlq5cSaV/y88Pw1HpVJhKddxFcaayJatHyqVig4pHWxxaS6gzLd1GhC0fCpar0mmOgrqWS3IJoP9eEH9MBqM9iVOEJ40gSPCSPMwFcaYvYJ3eUyaWO5sOA6etHJjTRhp9adTWiek/JoLataIAax6q1jqbCUEBQWhVKqQTLV70ir2i5IsrmE0Ge2GGYBaKYw0QSXCSPMwdo9ZuQcNq2hI3FhUKhVmq02xieVO10lLS0M2uZA8UJ40IIy01oEkSYSEhiKZSmo/rnx/eHi4J8RqMZhNZlRVPWnl4aEmU/26PAhaNsLd4GG0Wi0nTpwgfdvjKHVBtOnSlxjEcmdjUKlUWKy2XAyL8KS5TGpqKgBSvoQcVLc7TcoXSQOtjajICLIybFZ84Yk9nFn7LabifMI7DyZh+DVICoUw0hqIyWxGKeqkCZwgfsk8zMqVKzl58iQAJn0Bh9adJXzAAOH5aQRGo5ETJ05w6zd5dIoPQZaDhJHmAu3atUOpVGLNs0JSPU7Ih6joKJE00IqIjo7mwKls9GV6ji54F2v50mbWliVogsKJ7jsGyVhiP1ZQfywWC/tO5XL7P9nEB6v518BI+3aBQPySeZh//vnH4b1stVJSUmJrsitoEB988AEnT57kJLB8fxHJyUphpLmAWq0muW0yR/KPINcjME1VoCKtR5oHJBN4C1FRUUhGPSUZx+wGWgVFJ/eWG2l6FAqFiElzkbPp6fy+84j9/R+HionvHCV6dwoAEZPmcXr27Fltm3jybDi5ubns2LHDYVtWVpYw0lykY2pHlIX1KKhsBmuh1b5EKmgdREdHI5sN6MJjkRSO3xP/mLYASMZiwsIjRGFuF0k/m+7w/miOgaKiomaSRuBtCCPNw9x8883ExMQgSQqUukDa9xxIVJQw0hpKQEBAtYbOarVaGGkukpKSgrXUCnWVwiqsPF7QeoiJiQFAo1KQfMktqPyDQZIITbuAmP6XAqAw6ImLjWlOMX2S8+ORJWzZ/mJ1RQBiudPjhIWF0bFjR5JG/R+W2G74HVxGYKBoSNxQtFotU6dOZfbs2QDoNIrKAq2CetOhQwfbiwKglkTjiqQB+/GCVkGFkaYwFBPeZTBhnQciW8woVJWlg1QmPbGxXZtLRJ8lNTWFwrwsCkts2ZxTB0VyWOUnPJICQBhpHicgwGaQKZGxKhQoLEYCzvMECVxjwoQJbNmyhf9LziYgOIhPDoYIT5qLVBhdUr6EHFNLLEwBaP20xMbGekgygTdQ8feWDLZlOElSIFUx0LBakQ1FxMXFNYd4Pk1IcDAzxvWkW8A54oLVRAT5ceefooyQwIZY7vQwGo0GhVIJZluhQoXVRGBgYDNL5dtoNBq0Wi392gWhLjfORHFg1wgNDSU4JNi+nOkMqVCiXbt2KBRCdbQmIiIiUGs0KAw1x0pJxmKQZWGkNQCNRoMsKRjaPpAOkVp7H2JRlkkAwkjzOJIk4e8fgGSxGWmSxWj3rgkaRoVBZrZK9s4D4inUdTp06ICisHaVoCxS0r5dew9JJPAWJEkiNjbO7kmrtr98e0JCgifFahFo/fwwWCrjz4zllTfEg6YAhJHWLAQGVhppmI3Ck9ZIKpSZyYp4Cm0E7du1RyqqpT2UAaylVtq3F0ZaayQpMQGVsWYjTVFmc8HGx8d7UqQWgX9AAGVVjLTS8tfnJ0QJWifCSGsGggKDwGywxXFYxHJnY6kw0owWCZPVcZug/rRt29bWHspZH+3CyuMErY+EhASksiJba4/zkAyFqFQqoqKimkEy3yYgIJBSS+VPcalZKt8uVlgEwkhrFoKDg1BYTVDuTfv/9u43pqmz7wP49zotlWJpgWIpYHEx0eh8NjXL5nRTY42yxTENyCvjoi80Wfw74oxur6aLm4uLPs8LzTAxmhj3wjhlukw3XTYcW26zZWYhMRkvtig8CpkihUJ7+ufcL8o5FNHbe/TQ61C+nzfQP8CvBK5+z/WXO7dnJr0nTU3w7M7RmjZtWuqTJ8xLEyEx/Hk0oVRWVkJLxCBiI1O8EgnBX17BuYqj4HK50B8fWskZHgxpvHgngCFNCpfLBSWhQiSixm0aPaMnLSkQSwrk5dn5ZjEKgUDqTCjR+4T9mXqBPEcee0smqKlTpwIARKRnxGN2tReBqZyPNhputxt98aH/uXBMGPcT8Z1MApfLBZFQIQZXeDKkZSY/P7WxVzQhoCaBSRzqHJWSkhI4C5zAEzY7F70CgUCAAXiC0hcF6PPPDJoGEek1Qhz9M263G32qhuTgKHJfTDHuJ2JrK4HL5UrNSRsc7uTcg8zoiwRiyVRQy8/nooHREEKgKlAF0ff4njRbvw3TqjjUOVGVlZVBsdkgosNDmoj1Q0vEGNJGyePxQNOGetBCqkCBM5/zagkAQ5oULpcLWiJuzO3gnLTMpPekRRMCkxwMaaMVCARgCz9mp/MkkOxL8o14ArPb7Sgr84/oSRODt7n9xujoB9L3GiFNQVGRR2ZJZCEMaRLow5tKtG/YbRodvSdNTQioCSDf6ZRc0fg1depUJMNJIPnIA/0AknwjnuiqAlNhe2QbDj20McCPTklJCQCgR1UGPwoUl3hllkQWwpAmgR7KhMqQZgaHwwFFCEQTQDQpGNIyUFFRkdonLfzIA31pj9OEVVlZmQpladtwiGgIis0Gn88nsbLxS+9J61EHe9LidpQwpNEghjQJjJAW7YUQAk6GiowIkZqHlhruVOB0chPI0TJC2CMhTYTF8MdpQqqoqIAWV1NzagcpkRB8vjKelztKek9aKKr3pClGcCNiSJNAXyigqGE4Cwq4Ws4Ezvx8RBIC0aRizFGjf04/e1EPZYZwak5SaWmphKrIKowVnmmLB2xqH6oCHOocLbfbDZtNQbeaOtauN6oZwY2I6UCCoZ60MI/+MEn+4Pl30YTg7zQDXq831SPy6HBnGCj1lfKCYoIzQnxkaF6aEu3lweoZUBQFxUUe9EQVhAaHPL1eDndSCltcCfQQIZIxuCZzPpoZnAUFiCT0LTjYkzZaiqKgdEppaqFA+v39CirLuWhgotPDmKIftB5XocUiDGkZKvGWokdVjMUD7EkjHUOaBOn7onGPNHM4nQWIJgQiCXCOX4bK/eVQ+oc3DUpEQVlZmaSKyCqcTicK3R6IwZAmBld6+v1+mWWNeyUlXvTEbHg42JPGkEY6hjQJ0ofjXC6GNDOketIUROMaQ1qGfD4flEha05AEkv1Jrt4jAEC532+sTNe3EWJIy0xJSQlCMZuxeIALB0jHkCZBXl4ebLbUSigOzZkjPz8focHjVBjSMuPz+ZDsT6a24gCASOoDz+wkAPD7y2CPpcbDRTQ1eZG9rJkpLi5GT1QzhjsZ0kgnZc30w4cP8c4776CjowOVlZU4cuQIPJ7hOyxHo1GsW7cOqqoikUiguroa27dvl1HumJiUn4/+cB8nuZvE6XSiJ5oaKmDwzcyUKVNSAS0CwAljfhpXdhKQCmQi+nPqzE61DzabnaEiQ8XFxUgkgXv9CvInOXihSQYpPWmNjY1YuHAhvvnmGyxcuBCNjY0jnuNwOHDq1Cl8+eWXuHDhAq5fv46bN29mv9gxov8TMlCYw+l0IpYc+pxGz1hZNtiDxp40SjdlyhRoiRiQUKGoYZR4vVz1myG9k+Juvw1FHh4JRUOk/Gddu3YNa9asAQCsWbMGV69eHfEcIYQxqT4ejyMej0OIxx/8PB7p4YyBwhzpYZfBNzNGSEsdLQsR4bYANEQP60Lth4j1wzeFPayZ0kNaR1iBu6hIbjFkKVJC2v37941JyD6fDw8ePHjs8xKJBFavXo1FixZh0aJFmDt3bjbLHFP6eZP6R8oMQ5p59DCmhzMMpLbmeHRKAk1Mxt9HrB+2+AB7WE3gdrsBAP1xBR5PkdxiyFLGbE7ahg0b8Pfff4+4f+fOnf/197DZbGhqakIoFMKWLVvwxx9/YObMmSZWKY8y+P7HQGGO9B5J/k4zU6Rfyesn/0SBQk8hh7QIwND2ECI2AKEOsIfVBHpIe/RzojELaSdPnnziY16vF11dXfD5fOjq6nrqnjButxsLFizA9evXcyak6UO3DodDciW5Ib1HkiEtM5MmTYKzwIlwJLVyT0QESoq5bxOl6O21Eu2DFo9y0YAJ9FNoAKCwsFBiJWQ1Ui6Ng8EgLly4AAC4cOECli9fPuI5Dx48QCiUOh8uEongp59+wvTp07NZZlZwuNMcHO40l6fIY/SkCZUhjYZMnjwZNpsdykA3AG4XYYb0kMYNzimdlJC2efNmtLS0YOXKlWhpacHmzZsBAJ2dndi0aRMAoKurC2+99RZqamqwdu1aLFq0CMuWLZNR7phiT5o50sMuf6eZKykugRjc0kRRFb4Rk0EIgUK3G2LgIQBwrqIJ7Ha7MbrCkEbppOyTVlxcjFOnTo24v6ysDMePHwcAzJo1y+hty0WaltopNC8vT3IluSE9mLEnLXNFniIo/68giSQQ5RsxDVfk8aD79h0AnENlNu6dSek4E1gS/arJbpeSk3NOek8ah5Az53a7gQFA3BZIqkm+EdMwHo8HIhk3PqfM6Rfu3JaJ0jGkScahOXMwpJmrsrISWkSD8q9UE1FRUSG5IrISt3tocjsnupvjf+Y8C4DnoNJw7MaRpKqqCm1tbbxqMkn6sLHNZpNYSW5Yv349gsEgkskk7HY7QxoNkz5vinOozHHkf/8P/f39Q1vgEIEhTZodO3YgGAxizpw5skvJCek9krl0MoUsiqIgEAjILoMsSg9miqJwDqhJHA4HR1ZoBIY0SYqKirB48WLZZeQMLsAgyh49pE3Kz+dFEdEY4pw0ygkMaUTZo0/TsCmcWkA0lhjSKCdwmIAoe8rLywGAQ+JEY4zDnZQTuJUJUfYEg0E8//zzXNlJNMb4zkY5gSGNKLtKS0tll0CU8/jORjlh8uTJqKur4/ALERHlDIY0yglCCOzYsUN2GURERKbhwgEiIiIiC2JIIyIiIrIghjQiIiIiC2JIIyIiIrIghjQiIiIiC2JIIyIiIrIghjQiIiIiC2JIIyIiIrIghjQiIiIiC2JIIyIiIrKgnDwWqqOjA7W1tbLLIKIs6ejokF2Cadh+EU08T2rDhKZpWpZrISIiIqKn4HAnERERkQUxpBERERFZEEMaERERkQUxpBERERFZEEMaERERkQUxpEnQ3NyM6upqrFixAo2NjbLLGff27t2LhQsX4o033pBdSs64e/cu1q9fj9dffx2rVq3CqVOnZJdEFsI2zFxsw8yVU+2XRlkVj8e15cuXa7dv39ai0ahWU1OjtbW1yS5rXLtx44bW2tqqrVq1SnYpOaOzs1NrbW3VNE3Tent7tZUrV/LvlDRNYxs2FtiGmSuX2i/2pGXZ77//jmnTpiEQCMDhcGDVqlW4du2a7LLGtRdffBEej0d2GTnF5/Nhzpw5AACXy4Xp06ejs7NTclVkBWzDzMc2zFy51H4xpGVZZ2cn/H6/cbusrGzc/vHQxNDe3o5bt25h7ty5skshC2AbRuPJeG+/GNKyTHvMAQ9CCAmVED1dOBzG9u3b8d5778HlcskuhyyAbRiNF7nQfjGkZZnf78e9e/eM252dnfD5fBIrInq8WCyG7du3o6amBitXrpRdDlkE2zAaD3Kl/WJIy7LnnnsOf/31F+7cuQNVVfHVV18hGAzKLotoGE3T8P7772P69OnYuHGj7HLIQtiGkdXlUvvFA9Yl+OGHH3DgwAEkEgnU1dXh7bffll3SuNbQ0IAbN26gu7sbXq8X27ZtQ319veyyxrVffvkF69atw8yZM6EoqWu5hoYGLF26VHJlZAVsw8zFNsxcudR+MaQRERERWRCHO4mIiIgsiCGNiIiIyIIY0oiIiIgsiCGNiIiIyIIY0oiIiIgsiCGNLOfgwYN47bXXUFNTgy1btiAUCgEAWlpaUFtbi5qaGtTW1uLnn382vubw4cNYunQp5s+fP+x7qaqKnTt3YsWKFaivr0d7e/uwx/v6+rB48WLs27fPuG/Pnj0IBoNYvXo1Vq9ejVu3bo3hqyWiXMM2jMzCkEbSaJqGZDI54v5XXnkFly5dwsWLF/HMM8/gs88+AwAUFxfj2LFjuHjxIj7++GPs3r3b+Jply5bh7NmzI77X2bNn4Xa78e2332LDhg04dOjQsMePHDmCl156acTX7d69G01NTWhqasLs2bMzfalElIPYhtFYs8sugCaW9vZ2bNq0CQsWLMDNmzdRWFiI7u5uCCFQV1eHDRs24NVXXzWeP2/ePFy+fBkA8Oyzzxr3z5gxA6qqQlVVOBwOzJs377E/77vvvsPWrVsBANXV1di3bx80TYMQAq2trbh//z4WL16M1tbWsXvRRJQz2IZRNrEnjbLuzz//xJo1a/Dhhx/CbrcbV5y1tbUjnnvu3DksWbJkxP1XrlzB7Nmz4XA4/uPP6uzsRHl5OQDAbrcbDWoymcTBgweHXcmmO3z4MGpqanDgwAGoqjqKV0lEuYptGGULQxplXUVFBebNm4dAIIA7d+5g//79aG5uhsvlGva8Y8eOwWaz4c033xx2f1tbGw4dOjRsDsaTPO5ADSEEzpw5gyVLlhiNX7qGhgZcvnwZ586dQ09PDxobG//hKySiXMY2jLKFw52UdQUFBQAAj8eDpqYm/Pjjjzhz5gy+/vprfPTRRwCA8+fP4/vvv8fJkychhDC+9t69e9i6dSsOHjyIqqqqp/4sv9+Pu3fvwu/3Ix6Po7e3F0VFRfjtt9/w66+/4vPPP0c4HEYsFkNBQQF27doFn88HAHA4HKitrcWJEyfG4LdAROMV2zDKFoY0kubBgwdwOByorq5GVVUV9uzZAwBobm7G8ePHcfr0aTidTuP5oVAImzdvRkNDA1544YX/6mcEg0GcP38e8+fPx5UrV/Dyyy9DCIFPP/3UeM4XX3yB1tZW7Nq1CwDQ1dUFn88HTdNw9epVzJgxw8RXTUS5gm0YjTWGNJKmq6sLe/fuNVZHNTQ0AAD2798PVVWxceNGAMDcuXOxb98+nD59Grdv38bRo0dx9OhRAMCJEyfg9XrxySef4NKlSxgYGMCSJUtQX1+Pbdu2Ye3atXj33XexYsUKeDweHD58+Kl17dq1C93d3dA0DbNmzcIHH3wwRr8BIhrP2IbRWBPa4wa8iYiIiEgqLhwgIiIisiCGNCIiIiILYkgjIiIisiCGNCIiIiILYkgjIiIisiCGNCIiIiILYkgjIiIisiCGNCIiIiIL+jdi8RtMc1uFUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# genepair = 'RP1-29C18.10;ZNF501'\n",
    "# genepair = 'CCDC15;UNC5B'\n",
    "# genepair = 'GSTM3;RP1-29C18.10'\n",
    "# genepair = 'MMEL1;SARS2'\n",
    "genepair = 'AC005076.5;ARHGEF19'\n",
    "gene1, gene2 = genepair.split(';')\n",
    "\n",
    "if datasetname == 'ng':\n",
    "    ut_celltype = data_sc[data_sc.obs[dataset.celltype_id]==celltype]\n",
    "else:\n",
    "    ut_celltype = data_sc[(data_sc.obs[dataset.celltype_id]==celltype) &\n",
    "                          (data_sc.obs[dataset.timepoint_id_col]==dataset.chosen_condition['UT'])]\n",
    "\n",
    "ut_celltype_df = pd.DataFrame(data=ut_celltype.X.toarray(),\n",
    "                              columns=ut_celltype.var.index,\n",
    "                              index=ut_celltype.obs.index)\n",
    "selected_expression_df, ut_zscore_df, ut_zscore_p_df = get_individual_networks_selected_genepairs(\n",
    "    data_df = ut_celltype_df,\n",
    "    data_sc = ut_celltype,\n",
    "    individual_colname = dataset.individual_id_col,\n",
    "    genepair = genepair,\n",
    "    fillna=False\n",
    ")\n",
    "\n",
    "ut_t = ut_zscore_df.T\n",
    "ut_t['gt_sampleid'] = [gte_dict.get(name) for name in ut_t.index]\n",
    "ut_t = ut_t.set_index('gt_sampleid')\n",
    "common_individuals = list(set(gt.columns) & set(ut_t.index))\n",
    "gt_t = gt[common_individuals].T\n",
    "gt_t['genotype'] = [item.split(':')[0].count('1') for item in gt_t[0]]\n",
    "concat_df = pd.concat([gt_t, ut_t], axis=1).replace([np.inf, -np.inf], np.nan).dropna()\n",
    "print('Not Imputed', spearmanr(concat_df['genotype'], concat_df[genepair]))\n",
    "\n",
    "ut_t_imputed = ut_zscore_df.fillna(0).T\n",
    "ut_t_imputed['gt_sampleid'] = [gte_dict.get(name) for name in ut_t_imputed.index]\n",
    "ut_t_imputed = ut_t_imputed.set_index('gt_sampleid')\n",
    "common_individuals_imputed = list(set(gt.columns) & set(ut_t_imputed.index))\n",
    "gt_t_imputed = gt[common_individuals_imputed].T\n",
    "gt_t_imputed['genotype'] = [item.split(':')[0].count('1') for item in gt_t_imputed[0]]\n",
    "concat_imputed_df = pd.concat([gt_t_imputed, ut_t_imputed], axis=1).replace([np.inf, -np.inf], np.nan).dropna()\n",
    "print('Imputed', spearmanr(concat_imputed_df['genotype'], concat_imputed_df[genepair]))\n",
    "\n",
    "# dosage_dict = gt_t['genotype'].T.to_dict()\n",
    "# selected_expression_df_withsample = pd.concat([selected_expression_df,\n",
    "#                                                ut_celltype.obs[[dataset.individual_id_col]]],\n",
    "#                                              axis=1)\n",
    "# selected_expression_df_withsample['gt_sampleid'] = [gte_dict.get(name) for name in\n",
    "#                                                     selected_expression_df_withsample[dataset.individual_id_col]]\n",
    "# selected_expression_df_withsample['genotype'] = [dosage_dict.get(gt_sampleid) for gt_sampleid in\n",
    "#                                                selected_expression_df_withsample['gt_sampleid']]\n",
    "\n",
    "sns.set_style('white')\n",
    "refallele = gt['REF'].values[0]\n",
    "altallele = gt['ALT'].values[0]\n",
    "snp_name = f'{snp_id}_{altallele}'\n",
    "\n",
    "_, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\n",
    "ax1, ax2 = axes\n",
    "\n",
    "im_coef, im_p = spearmanr(concat_imputed_df['genotype'], concat_imputed_df[genepair])\n",
    "sns.violinplot(x=concat_imputed_df['genotype'],  \n",
    "               y=concat_imputed_df[genepair], \n",
    "               ax=ax1,\n",
    "              inner=None)\n",
    "sns.swarmplot(x=concat_imputed_df['genotype'],  \n",
    "               y=concat_imputed_df[genepair], \n",
    "               ax=ax1,\n",
    "              color='black')\n",
    "ax1.set_title(f'Imputed r={im_coef:.2f}; pvalue {im_p:.4f}')\n",
    "# ax1.set_xticklabels([f'{refallele}{refallele}', \n",
    "#                      f'{refallele}{altallele}',\n",
    "#                      f'{altallele}{altallele}'])\n",
    "ax1.set_xlabel(snp_id)\n",
    "\n",
    "coef, p = spearmanr(concat_df['genotype'], concat_df[genepair])\n",
    "sns.violinplot(x=concat_df['genotype'],  \n",
    "               y=concat_df[genepair], \n",
    "               ax=ax2,\n",
    "              inner=None)\n",
    "sns.swarmplot(x=concat_df['genotype'],  \n",
    "               y=concat_df[genepair], \n",
    "               ax=ax2,\n",
    "              color='black')\n",
    "ax2.set_xlabel('')\n",
    "ax2.set_title(f'Not Imputed r={coef:.2f}; pvalue {p:.4f}')\n",
    "# ax2.set_xticklabels([f'{refallele}{refallele}', \n",
    "#                      f'{refallele}{altallele}',\n",
    "#                      f'{altallele}{altallele}'])\n",
    "ax2.set_xlabel(snp_id)\n",
    "plt.savefig(example_savedir/f'{snp_name}_ref{refallele}_alt{altallele}_{gene1}_{gene2}.{celltype}_{datasetname}.full.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/tools/Beeline/miniconda/envs/scpy3.8/lib/python3.8/site-packages/seaborn/categorical.py:1296: UserWarning: 42.5% of the points cannot be placed; you may want to decrease the size of the markers or use stripplot.\n",
      "  warnings.warn(msg, UserWarning)\n",
      "/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/tools/Beeline/miniconda/envs/scpy3.8/lib/python3.8/site-packages/seaborn/categorical.py:1296: UserWarning: 7.1% of the points cannot be placed; you may want to decrease the size of the markers or use stripplot.\n",
      "  warnings.warn(msg, UserWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAFNCAYAAACAH1JNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABur0lEQVR4nO3dd1zV9ffA8dcdXLaAIA5wY2ruvbUwNQc5UtNsqBmapqmVXyuzMrMsy/qaZVjatBwpZurXnXvk+pkzcaSiguzNXZ/fH1duXhmiApd7Oc/HA+Uz77kX7ptz31OlKIqCEEIIIYRwGGp7ByCEEEIIIe6OJHBCCCGEEA5GEjghhBBCCAcjCZwQQgghhIORBE4IIYQQwsFIAieEEEII4WAkgRN3VLduXf755x97h+EUpk6dyty5c+0dhhAiH/v376dz5872DsNphIaGsmfPHnuH4ZQkgSsh9volXrlyJUOHDi3xx71Xa9as4eGHH6Zp06aMHTuWpKSkPM+Lj49n8uTJdOzYkRYtWjBkyBD+7//+z3p8//791KtXj2bNmlm/Vq1aVULPwj4K+9oBXLlyhaeffpomTZrw6KOP2vxu7tu3j7CwMFq2bEmbNm0YN24cMTEx1uNJSUlMnDiRNm3a0KZNG15++WXS0tIASEhIYMiQIbRp04aWLVvyxBNPcOjQoWJ7zuLehIaG0r59ezIyMqz7li9fztNPP12o659++mmWL1+e7/ErV65Qt25djEbjfcd6txzpQ5KiKHz00UfW99KHH35IQVOz7t27l0cffZQmTZrw9NNPEx0dbT327bff0rVrV5o3b07Hjh2ZNWuWXV7/knK3r926devo2bMnzZo1o1evXmzevNl6LCUlhf/85z+0a9eOdu3aMW/evDzvceDAAerWrWvz+/XHH38wdOhQWrZsSYcOHZg2bZq1PCxuksCJIqMoCmaz+Z6vP3v2LNOnT+fDDz9k9+7duLu788477+R5bkZGBo0aNWLlypUcOHCA/v37Ex4eTnp6uvWcwMBAjhw5Yv3q37//PcdW2t3Nawfw8ssv8+CDD7J//34mTZrEhAkTSEhIACAkJISvv/6agwcPsnPnTqpXr85bb71lvfbTTz8lJSWFLVu2sHnzZuLj460FnqenJ7NmzWLv3r38+eefPP/887zwwgtO/YfEUZlMJr7//nt7h+HQ7vf3eunSpWzevJnVq1fz22+/8ccff/DLL7/keW5CQgIvvvgiL730EgcOHKBhw4ZMmjTJejw0NJRVq1Zx+PBhfv/9d06fPs0PP/xwX/GVZnfz2sXExDBlyhSmTp3K4cOHmTJlCi+//DLx8fEAvP/++2RmZrJ161aWL1/O6tWr+fXXX23uYTAYeO+992jSpInN/tTUVF544QV27tzJunXruH79Oh9++GHxPOnbSAJnBytXrmTIkCHMmjWLli1b0rVrVw4fPszKlSvp0qUL7dq1s6ktmjp1KtOnT2fEiBE0a9aMp556yvrJK69Pujmfjs+dO8dbb73F0aNHadasGS1btgRAr9cze/ZsHnroIdq3b8/06dPJysqyXv/111/TsWNHOnbsyIoVKwp8Lk8//TRz585lyJAhNGnShMuXL9/z67JmzRpCQ0Np1aoVnp6evPTSS2zatCnPTzNVq1ZlxIgRBAYGotFoeOKJJzAYDFy4cKFQjxUREcHo0aPzPR4aGspXX31Fr169aNWqFa+99hrZ2dkA9OzZk23btlnPNRqNtGnThhMnTgAwYcIEOnToQIsWLRg2bBhnz57N8zHyqh29tbn6Tj+nW93Na3fhwgVOnDjB+PHjcXNzo0ePHjzwwANs2LABgICAACpWrGg9X6PRcOnSJev2lStX6Nq1K15eXnh7e9OtWzeioqIAcHV1pVatWqjVahRFQa1Wk5ycTHJycr6vtbCP5557jkWLFpGSkpLn8cOHD/P444/TokULHn/8cQ4fPgzA3LlzOXjwIDNmzKBZs2bMmDHjjo81depU3n77bUaNGkWzZs0YMmQIN27c4L333qNVq1Y8+uijnDx50np+Qe+/gt43S5cuZc2aNXzzzTc0a9aMMWPGAJY/4OPHj6dt27aEhobaJK5ZWVlMnTqVVq1a0atXL/76668Cn0vdunX56aef6N69O927d7/jcy9IZGQkI0eOpFKlSlSsWJERI0bk21KwadMm6tSpQ8+ePXF1dWX8+PGcPn2ac+fOAVCtWjXKlSsHYH3v3dr1ZfTo0UREROR575y/I0uXLrWW/YsWLQIsr13jxo1tavRPnjxJmzZtMBgMXLp0iWeeecamRj6/36nba0dvb64u6Od0P6/d9evX8fb2pkuXLqhUKh566CHc3d2t5drWrVsZNWoU7u7uBAcHM3DgwFwJ3OLFi+nQoQO1atWy2R8WFkbnzp1xd3fHx8eHwYMHc+TIkXzjLkqSwNnJsWPHqFu3Lvv376dPnz5MnjyZv/76i02bNvHRRx8xY8YMm9qkNWvWMHbsWGvT4CuvvHLHx6hduzbvvPMOTZs25ciRIxw8eBCAjz76iAsXLhAZGcnGjRuJjY1l/vz5AOzYsYNFixaxaNEiNm7cyN69e+/4OKtXr+bdd9/l8OHDVKlShbfffpuWLVvm+RUWFpbvfc6ePUvdunWt29WqVcPFxYWLFy/eMYZTp05hMBioXr26dV9CQgLt27cnNDSUWbNm2TQXhYeH89VXXxV4z5w/BJs2beLChQt88cUXAPTu3Zvff//det6uXbvw8/OjQYMGAHTu3JkNGzawd+9eHnzwwUL9rPJS0M/pdnfz2kVFRVG1alW8vLys++rVq2dNwgCuXr1Ky5Ytady4MYsWLWLUqFHWY8OGDeOPP/6wJmYbNmygU6dONo8RFhZG48aNeeGFFxg0aBD+/v739BqI4tOwYUNat27NN998k+tYUlISo0eP5umnn2b//v2MGDGC0aNHk5iYyKRJk2jZsiXTp0/nyJEjTJ8+vVCPt379eiZOnMi+ffvQ6XQ88cQTNGjQgH379tGjRw/ef/99m/Pze/8V5IknniAsLIznnnuOI0eOsGDBAsxmMy+88AJ169Zlx44dfPfdd3z33Xfs3LkTgM8//5xLly6xadMmvvnmGyIjI+/4OJs3b2bZsmWsW7cOwNrlIK+vt99+O9/7nD17lnr16lm369Wrl+8Hvtvf4x4eHlSrVs3mfbtmzRqaN29O27ZtOX36NEOGDLEe++qrrwgPDy/wee3fv5+NGzfyzTffEBERwZ49e6hYsSJNmzZl48aNNo/To0cPXFxcUBSF0aNHs3PnTtavX8/169fzbYIsyJ1+Tnm9HoV97Ro2bEjt2rXZsmULJpOJzZs3o9PpbF7PWymKYnOv6Ohofv31V8aNG3fH5/Hnn38SEhJyx/OKgiRwdhIcHMzjjz+ORqOhV69eXLt2jXHjxqHT6ejYsSM6nc6m1uOhhx6iVatW6HQ6Jk2axNGjR7l27dpdP66iKCxfvpzXX38dX19fvLy8GD16NGvXrgUsheyAAQN44IEH8PDw4MUXX7zjPfv370+dOnXQarW4uLjw9ttvc/DgwTy/1qxZk+99MjIy8Pb2ttnn5eVlk8jmJS0tjSlTpvDiiy9ar69VqxaRkZHs2rWL7777jhMnTvDBBx/c8bncatiwYVSuXBlfX19eeOEF62sUFhbG1q1byczMBCyFWZ8+fazXDRw4EC8vL3Q6nfVTcmpq6l099p1+Tre7m9cuPT0917ne3t4251apUoWDBw+yb98+XnrpJZtPnQ8++CAGg8H6iVuj0fDkk0/a3G/NmjUcOnSIjz/+mBYtWtzVcxclZ8KECfz444/W5vMcf/zxB9WrV6dfv35otVr69OlDrVq1bGqe71a3bt1o2LAhrq6udOvWDVdXV/r162ctA0+dOmVzfn7vv7v1119/WZsfdTodVatWZfDgwdbka/369YwZMwZfX18qV65cqH6A4eHh+Pr64ubmBlh+3/Mr8wpK4DIyMmw+SHl7e5ORkZFnX67CvMfDwsI4fPgwGzZsYMiQIXf9wWncuHF4eHhQt25dBgwYYP2gGhYWZv1eURTWrVtn/TBevXp1OnTogE6no3z58owYMYI///zzrh4X7vxzut3dvHYajYa+ffvyyiuv0KhRI15++WVmzJiBh4cHAJ06dSIiIoK0tDT++ecffv31V2v5DjBz5kxeeuklPD09C3wOu3fvJjIykgkTJtz1878X2hJ5FJHLrW+snEIgICDAus/V1dXmjVmpUiXr956envj4+BAbG3vXb9CEhAQyMzMZMGCAdd+tfddiY2Np2LCh9VhQUNAd71m5cuW7igHg4MGDPP/884AlWVi7di0eHh65mvzS0tIKfNNkZWUxZswYmjRpYtMkWqFCBSpUqABYmltfffVVRo8eXajmnhy3Pq8qVaoQGxsLWAqs2rVrs23bNh5++GG2bt1q/dRuMpmYO3cu//vf/0hISECttnxGSkxMzFX4FuROP6fb3c1r5+npWehzfX196d+/P3379mXHjh1otVpeeukl6tWrxxdffIGiKMyePZtXX32Vzz77zOZaV1dX+vTpQ8+ePalfv77Np2VROjzwwAM89NBDREREULt2bev+2NhYqlSpYnNulSpVbAaz3K3by7xbyzs3NzebGnLI//13t6Kjo4mNjbV2IQHL+zRnOzY2Ntdj3cm9lHkLFiyw1vqHhYVZE4hby/m0tDQ8PDxQqVS5rs/rPZ6enp7n+7ZGjRrUqVOHd955h88//7zQMd76vIKCgvj7778B6NGjB++++y4xMTH8888/qFQq6+sXHx/PzJkzOXjwIOnp6SiKYm3KvRt3+jnd7m5euz179jBnzhy+//57GjRowPHjxxk7diwLFy6kfv36TJs2jXfffZcePXrg6+tL7969rR8Ytm7dSnp6Or169Sow/qNHj/Lyyy/z3//+l5o1a971878XksA5iOvXr1u/T09PJzk5mcDAQFxdXQFLIpPzaeTGjRvWc2//Zfbz88PNzY21a9fa9HPKERgYaFOzd/Xq1TvGdvtjTJ8+Pd+atpxkrWXLlrn6CdSpU4fTp09bty9fvozBYKBGjRp53kuv1zNu3DgqVqx4x8RMpVIVOEIpL7e/DoGBgdbtPn368Pvvv2M2mwkJCbE23a5Zs4YtW7awePFigoODSU1NpVWrVnk+tru7u02ftlt/bnf6Od3ubl67kJAQLl++TFpamvV35vTp0za1iLcymUzEx8eTlpaGr68vZ86c4e2337Z+eh06dGiuGrhbGY1GLl++LAlcKTVhwgT69+/PyJEjrfsCAwNzvfevXbuWq6m8OOX3/ivofQO5y6PKlSsTHBxs0wR4qwoVKnDt2jXq1KmT63Hzc/tj9O7dO9+yMidZGzNmjLVPXo6c923jxo0By/swJ47b1alTx6aPV0ZGBpcuXcq3uc5oNNq04hTGtWvXrIn8ra95uXLl6NChA+vXr+f8+fP07t3b+hp8/PHHqFQqfvvtN/z8/Ni8eXO+5fHtP7u4uDjr93f6Od3ubl67U6dO0bJlSxo1agRA48aNady4MXv27KF+/fr4+vry8ccfW8//5JNPrPfdu3cvx48fp0OHDoBl0IJGo+Hvv//myy+/BCx9Al944QVmzZpFu3btChV/UZAmVAexfft2Dh48iF6v57PPPqNJkyZUrlyZ8uXLU7FiRVavXo3JZGLFihU2Awn8/f2JiYlBr9cDoFarGTRoELNmzbKOwImJibH2M3j00UdZtWoVUVFRZGZm3tWntxwzZsywGf1561dBzSBhYWFs27aNgwcPkpGRwWeffUa3bt1sqslzGAwGJkyYgKurK7Nnz7bWdOXYv38/V69eRVEUrl27xpw5c+jatav1+Lx58+7YVLJkyRKuX79OUlKStUN1jl69erF7925+/vlnm8QnPT0dnU6Hn58fmZmZfPLJJ/neP6fPxqlTp8jOzrbpN3Knn9P9vHY1a9akfv36zJ8/n+zsbDZt2sSZM2fo0aMHABs3buT8+fOYzWYSEhJ4//33efDBB/H19QUs/UmWL19OVlYWWVlZLF261NqX5OjRo9bf06ysLCIiIoiLi7MWhqL0qV69Or169bIZsdilSxcuXrzImjVrMBqNrFu3jqioKB566CHA0lpwPwOWCiO/919B7xuwlHlXrlyxbjdu3BgvLy8iIiLIysrCZDLx999/c+zYMcAyKCkiIoLk5GSuX79+TyM3165dm2+ZV9CHy759+7J48WJiYmKIiYlh8eLF+Y6W79atG2fPnmXDhg1kZ2czf/586tata024li9fbi0roqKiiIiIsEkmnn766Tv2Tfviiy/IzMzk7NmzrFy50qbMCwsLY/Xq1WzYsMGmL3N6ejoeHh6UK1eOmJgYvv7663zvX79+fbZv305SUhI3btzgu+++sx6708/pfl67Ro0acfDgQWsz/cmTJzl06JC13Lp06RKJiYmYTCa2b9/O0qVLeeGFFwB46aWX2LBhA5GRkURGRhIaGsqgQYOsfTb//vtvRo0axZtvvkloaGiBr29RkwTOQfTp04f58+dbRzt+9NFH1mPvvvsu33zzDW3atCEqKopmzZpZj7Vt25aQkBA6duxImzZtAHj11VepXr06gwcPpnnz5gwfPtw6erNLly48++yzPPvss3Tr1o22bduW2HPMqfJ/5ZVXaN++Penp6TbTV0yfPt3aYfrIkSNs27aN3bt306pVK+tcbzkDNU6ePMkTTzxB06ZNGTJkCA888ABvvPGG9V7Xrl2jefPmBcbTp08fRo4cySOPPELVqlWtb2iw1FDkDA65tZDr168fVapUoVOnTvTu3ZumTZvme/+aNWsybtw4hg8fTvfu3XP1FSvo53Q/rx1YPmEeP36cVq1aMWfOHP773/9Svnx5wJIojho1iubNmxMWFoZarbZJ5GfNmkV0dDRdunShc+fOXL582dq/UK/XM2PGDNq0aUPnzp3ZsWMHERERhapFFPYzbtw4myZMPz8/FixYwOLFi2nTpg1ff/01CxYssP6OPPPMM2zYsIFWrVoxc+bMYokpv/ffnd43AwcOJCoqipYtWzJ27Fg0Gg1ffvklp0+fpmvXrrRt29Zmrq4XX3yRKlWq0LVrV0aOHEnfvn2L5fnkZciQITz88MOEhYURFhZGly5dbAYe9O7dm99++w2A8uXLM2/ePObOnUurVq04duyYzQfEw4cPExYWRtOmTQkPD6dz585MnjzZerwwZV7r1q3p1q0bw4cPZ+TIkXTs2NF6LDQ0lIsXLxIQEGBTm/7iiy9y8uRJWrZsSXh4eIEjc/v27Uu9evUIDQ1l5MiRNmXnnX5O9/PatW7dmvHjxzNhwgSaNWvG+PHjGT16tPX5HT9+nLCwMJo3b84nn3zCnDlzrLV5Xl5e1i45FSpUwM3NDXd3d+sH2sWLF5OQkMAbb7xh/TvUu3fvAl/noqJS7rZdSZS4qVOnUrFiRZs5f8T96du3L99++y1+fn55Hg8NDWXmzJm0b9++hCMTQsj7r2hdv36dl156iaVLl+Z5PGdqoBMnTqDVSs8qRyE/KVEmrV692t4hCCFEiahUqVK+yZtwXNKEKoQQQgjhYKQJVQghhBDCwUgNnBBCCCGEg5EETgghhBDCwZSpQQxt2rQp1MoCQgjnER0dzf79++0dxn2T8kuIsqeg8qtMJXBBQUGsXLnS3mEIIUrQrcuROTIpv4Qoewoqv6QJVQghhBDCwUgCJ4QQQgjhYCSBE0IIIYRwMJLACSGEEEI4GEnghBBCCCEcjCRwQgghhBAORhI4IYQQQggHY9cEbseOHfTo0YNu3boRERGR67iiKMycOZNu3boRFhbGiRMnCn3t3frzzz/p1asXmZmZXL58mYyMDJvjN27cID4+3mZfWloaV65csdlnNBq5cOECJpPJZr+j3XPOnDlMnz6d+xEbG0tCQgIZGRlcvnzZ5pjJZOLChQsYjUab/VeuXCEtLc1mX3x8PHFxcTb7HO2eiqJw8eJF9Ho9ji4zM5NLly7ZOwwhhCjbFDsxGo1K165dlUuXLinZ2dlKWFiYcvbsWZtz/vjjD+W5555TzGazcuTIEWXgwIGFvjYv/fv3z3O/l5eXAth8eXt7K4sWLVKMRqMyfPhwRa1WKxqNRhk7dqxiNpuVefPmKR4eHgqgdO7cWYmPj1d27typVKlSRQGUGjVqKAcPHlSuXbumtG7d2qHuefvXhQsX7upnazAYlGHDhikqlUrRaDSKi4uLAiht27ZVYmJilD///FOpVq2aAihBQUHK7t27lfj4eKVTp04KoHh6eirz589XzGazMmbMGEWj0ShqtVoZMWKEYjQalW+++Ubx9vZ2qHu6ubkpgBIYGKhs2LDhrl7P0mTJkiWKj4+PAijNmzdXLl++bO+Q7ii/972jcZbnIYQovILe93ZL4A4fPqyMHDnSur1gwQJlwYIFNue8+eabypo1a6zb3bt3V2JiYgp1bV7yeiGuXLmSZ9KS80d34cKFufZ///33ikajsdk3ZcoUpW7dujb7WrdurYwePdoh73nrl6ur6139bL///vt87zV27FilRYsWNvsefPBB5eWXX7bZp9Vqle+++y7X9d98843i6urqkPfM+QoODlaMRuNdvaalQUpKSq4PO88++6y9w7ojZ0l8nOV5CCEKr6D3vd2W0oqJiaFSpUrW7YoVK3Ls2LECz6lUqRIxMTGFurawvvnmm3yPZWVl5bkG2Z49e3I1PZ44cYK///7bZt/Jkydxc3MrkXseOHCgSO95q+zs7HyP5eXkyZP5Hjt16hSnTp2y2Xf69GmqV69us89oNLJ3795c1+/fvz9XPCV5zxo1atzzPXNcuXKF1NRUfH198zxeWl2+fDlXs3FBP2shhBDFx2594BRFybVPpVIV6pzCXFtYBfXzqlq1KsOHD7e5t0aj4bnnniMgIMDm3LCwMB599FGbfX369KF3794lcs9nn322SO95q9sToTsp6F69e/fOdbxXr1706dPHZl+FChV47rnnUKv//RVVqVSMGDEi14LeJXnP2/fdzT1zdOjQweGSN4B69epRu3Ztm323v8ZCCCFKSElVA96utDShKoqiPPnkkzbNQkFBQcpjjz2mHD9+XFEURVm6dKnSpk0bpUOHDsrq1asVRVGUQ4cOKT179lQaN26szJ49WzGbzcqNGzeU4cOHK/Xr11dGjx6tJCUlKUajUXnnnXeURo0aOcw91Wr1PTef5vjpp5+U1q1bK82bN1dat26tNGrUSHn33XcVk8mkJCYmKuHh4Ur9+vWVESNGKHFxcYrZbFbef/99pXHjxkqvXr2UI0eOKIqiKJGRkUr79u2Vtm3bKsuWLVMURVH++usvJSwszGHu2aBBA6Vt27bKgw8+qDz55JNKdHT0Pb2mpcGZM2eU/v37Kw0aNFCmTZumGAwGe4d0R87S9Ogsz0MIUXilsg+cwWBQQkNDbQYi/P333zbnbNu2zWYQw+OPP17oa/MiBaAQZY+zvO+d5XkIIQqvVPaB02q1TJ8+nVGjRmEymXj88cepU6cOP//8MwBDhw6lS5cubN++nW7duuHu7s6sWbMKvFYIIYQQoiywWwIH0KVLF7p06WKzb+jQodbvVSoVb731VqGvFUIIIYQoC2QlBiGEEEIIByMJnBBCCCGEg5EETgghhBDCwUgCJ4QQQgjhYCSBE0IIIYRwMJLACSGEEEI4GEnghBBCCCEcjCRwQgghhBAORhI4IYQQQggHIwmcEEIIIYSDkQSuFMrSm4hPy8ZkVuwdihBCCCFKIbuuhSrypqCQkmkgPduIj7sL3m4uqNUqe4clhBBCiFJCauBKMaNZIT5dz7XkTDL1RnuHI4QQQohSQhI4B5BtNHM9OYu4tGzM0qwqhBBClHnShOogFCAl00Cm3oinqxYPFy1uOo29wxJCCCGEHUgC52AMJoWkDAPJGHDRqCnnrsXLVfrICSGEEGWJJHAOSgH0JjNxaXqSMgx4umrx1GlxdVGjUkkyJ4QQQjgzSeCcgNGskJxpICXTgEatwtNVi5ebFletNLEKIYQQzkgSOCeicEsyl2WgnJsLfh46aV4VQgghnIwkcE5KUSA500CG3oivhw4vV600rQohhBBOQhI4J2cwKdxIzSYpQ28ZvarT4uYiTatCCCGEI5MEroywjl7NMKDTqnHXaXDTanB10aCRJlYhhBDCoUgCV8YoWCYGzjaaUWFApQIXjRpXFw0eLhrcdRppahVCCCFKOUngyjAFS1+5nIQuNdMyt5yPh4v0mRNCCCFKMUnghFXO3HI3UrNJTNfj4arFVaNGq1Hh5iI1c0IIIURpIQmcyJPRrJCSaQBABWjUKjx0WrQaFTqNGjcXjUxPIoQQQtiJJHClzMrDV9h44jp1K3nTvnYAFbxd7R2SdX65lKx/EzpL3zkNri5qXDVqXDRqdFq1JHVCCCFECZAErpRZtPsCx6NT+N+JGD7bEkW9St50DAmgY0gA1fw97B0ecGvfORPZRhNgm9S5uVhq6Fy1arQatV1jFUIIIZyRJHClzCeDm/L51ii2nYklNcvI6eupnL6eyte7LlDVz52OdSzJXN1K3qhLUZ+025O65EwDahWoVZb+c+6uGty1GknohBBCiCIgCVwp80BFb2Y/3ogrCRn8X3Qyu87GsTsqnhtp2VxOzOTnA5f5+cBl/L10dKgdQMcQf5pU9cWlFCZGZgXMikJatpG0bCPq25tdtWp0Gml2FUIIIe6WJHCllFajpnk1P5pX82N8aAhnYlLZHRXPrqg4/onPID5Nz2//d5Xf/u8qnq4a2tXyp2NIAK1qlMddVzpXWjAX0OzqoVPjerOGzkWjkhGvQgghRAEkgXMAKpWKepXKUa9SOZ7rWJPLCRnsjopjV1QcJ6+lkp5tYvOpWDafisVFo6JFdT86hQTQrrY/vh46e4efr9x96SwTC6uwTC6s1ajRqlWWr5sDJSS5E0IIISSBc0hVy3swpHU1hrSuRlxaNnvPWWrmjlxKwmBS2Hc+gX3nE1CroGGQDx1CAugUEkAlHzd7h35HimK7WsStVDf71OUkdJqcbA/rf6hu2VarVKg1KlxUluZaWTJMCCGEs5AEzsEFeLkS1qQKYU2qkJZtZP/5eHZFxXPgQgKZBhPHriRz7EoyX/5xjpAKXnQI8adjnQBqBXg6XE2WooBJUTCZlVzJXUFymmq1astUJy4aNVqtChWWZBGz5b7KzT57Of9bvkC5+b1apUKtVuGiVqNWY0kgb4/xtnhv3avcdq5apUKjUaG9ed+c829NVO39M1Juvt5Gk4LBbAYFFBTUahVeri52jU0IIcoySeCciJerlq71K9K1fkX0RjOHLyWyKyqOPVHxJGUaiLqRRtSNNL7b+w+VfdzoGBJAhxB/GlTxceraqZymWr3JjN5U+MQv7ztBJqYiiStHnq+86t9aRI3akuBp1ap/t1X/JohqlQqVGtSoLM8VBbPZkniZzJZk1FJZefOmOYkq/yasys2bqVSWe6tUKgxGM0az2fKsFdsE1M1FIwmcEELYkSRwTkqnVdO2lj9ta/ljekTh5NUUdt3sN3ctOYtryVksP3SF5Yeu4OfhQrvalkEQzav5odOWvhGtzuz2mrmcnTkJltGc5xlCCCHKMEngygCNWkWjYB8aBfswpkstzselW6cnibqRRmKGgXV/XWfdX9dxd9HQumZ5Oob406aWP16u8isihBBClDby17mMUalU1K7gRe0KXjzbvgbXkjPZFRXP7qg4jkcnk2kwsf3vG2z/+wZatYqmVX3pWCeA9rX9CfCy/7JeQgghhJAErsyr7OPOoBbBDGoRTFKG/uaI1ngO/pOAwaRw8J9EDv6TyKebz/JgZW863FzWq2r50rGslxBCCFEW2SWBS0pKYtKkSURHRxMUFMSnn36Kj49PrvN27NjBe++9h9lsZtCgQYSHhwMwb948li1bRvny5QGYPHkyXbp0KdHn4Ix8PXT0bFSZno0qk6k3ceBiAruj4th3PoG0bCMnr6Vy8loqC3deoLq/h3WN1gcqetl9tKQQQghRltglgYuIiKBdu3aEh4cTERFBREQEr776qs05JpOJGTNmsHjxYipWrMjAgQMJDQ0lJCQEgOHDh/Pcc8/ZI/wywV2nocsDFejyQAWMJjNHLydZVoI4F0d8mp5/4jP4J/4SP+2/RAUvV+v0JI2DfGS9UyGEEKKY2SWB27JlCz/88AMA/fr14+mnn86VwB07dozq1atTtWpVAHr37s2WLVusCZwoOVqNmpY1ytOyRnnGdw3hzPVUy4jWs3FcTszkRlo2kUevEnn0Kt5uWtrV8qdDSAAta/jh7lI6l/USQgghHJldErj4+HgCAwMBCAwMJCEhIdc5MTExVKpUybpdsWJFjh07Zt3+6aefiIyMpGHDhkydOjXPJlhR9NQqFfUrl6N+5XI836kWl+IzrNOTnL6eSmqWkY0nY9h4MgZXrZqW1f3oWCeAtrX88XGXecOEEEKIolBsCdzw4cOJi4vLtX/ixImFul5Rcs99ldPPaujQoYwdOxaVSsVnn33GBx98wPvvv39f8Yp7U83fgyf9q/Fkm2rcSM1md1Qcu6PiOHolmWyjmd3n4tl9Lh61ChoH+9AxJID2IQFUKlf6l/USQgghSqtiS+C+/fbbfI/5+/sTGxtLYGAgsbGx1sEIt6pUqRLXr1+3bsfExFhr7QICAqz7Bw0axJgxY4oucHHPKni70q9ZEP2aBZGaZWDfecsgiAMXEsgymjl6OZmjl5P5fNs5QgK96BQSQMc6AdTw95BBEEIIIcRdsEsTamhoKJGRkYSHhxMZGUnXrl1zndOoUSMuXrzI5cuXqVixImvXruXjjz8GsCZ/AJs3b6ZOnTolGr+4M283F7o9WJFuD1Yk22Di4D+J7I6KZ8+5OFKyjETFphEVm8biPRcJ8nW3DIIICeDBKuVQSzInhBBCFMguCVx4eDgTJ05kxYoVVK5cmc8++wyw1LJNmzaNhQsXotVqmT59OqNGjcJkMvH4449bE7WPPvqI06dPAxAUFMSMGTPs8TREIbm6aOgQEkCHkABM5gc4Hp1s7TcXk5JNdFImyw5eYdlBy7JeOXPNNa3qK8t6CSGEEHlQKXl1NnNSAwYMYOXKlfYO444y9UauJ2flvUamE1EUhajYtJvJXDwX4tJtjnvoNLSpWZ6OIQG0rlkeT1nWq9Rwc9FQxdfd3mEUiqO87+/EWZ6HEKLwCnrfy19EYTcqlYo6Fb2pU9GbER1qEp2Uye6b05OcuJpCht7EtjM32HbmBi4aFc2q+dExxJ/2tQMo76mzd/hCCCGE3UgCJ0qNIF93BresyuCWVUlI17PnnGWN1sOXEjGYFA5cSODAhQTmbjrLg1XKWVaCqBNAkIPUBAkhhBBFRRK4UkilUuGu02JWFMxmBbOioPDv1CrKzX+cuYm1vKeOPo0r06dxZTL0Rg5cSGDn2Tj2X0ggQ2/ixNUUTlxN4asd56kZ4EnHEMvkwXUCZVkvIYQQzk8SuFLIzUVDJZ/cKxgoioJZAdPNpM5kVjCYzGQbzRhNZoxmBUVRyOnV6CwJnodOy0N1A3mobiB6481lvc7FsTsqnoR0PRfi0rkQl84P+y4R6O1qrZlrFOSDRi3JnBBCCOcjCZwDUalUaFQUmJSYzApGsxmzWcGoKJhMCtlGM9kGk6Umz8GTO51WTeua5WldszwvdVU4dS2FXWfj2H0uniuJmcSmZrPySDQrj0RTzk1Lu9qW6UlaVvfDVZb1EkII4SQkgXMyGrUKjTrvRCUnuTOZFFKyjGTqjQ6byIFlWa8GVXxoUMWH8M61uBifYRkEERXH3zFppGQZ2XAihg0nYnDTqmlVszwdQgJoV6s83m6yrJcQQgjHJQlcGWJN7rTg4aolPdtIapaRLIORWyeTccSkTqVSUTPAk5oBnjzVtjoxKVnsjopnV1Qcx64kkWU0s/NsHDvPxqFRq2hyc1mvDiEBVPB2tXf4QgghxF2RBK4M83TV4umqtfSfM+UMkLjZ5HqzX53BZMbsgBldxXJuDGgexIDmQSRnGth33pLMHbyYSLbRzOFLSRy+lMR/t0ZRt5I3HW+uBFHd39PeoQshhBB3JAmcQKtRo72l1dX9linWjCYzmXoT6XoTWQajQyZzPu4u9GhQiR4NKpFlMHHwYiK7z8Wx51w8qVlGzlxP5cz1VL7ZdZFgP3fLIIiQAOpV9pZlvYQQQpRKksCJAmk1arzd1Xi7u1iSOYOJbIMZvclSQ2e+OTLWUbi5aOhYxzJK1Wgycyw6mV1nLclcbGo2VxIz+eXPy/zy52X8PXW0v1kz17SqLy4aWdZLCCFE6SAJnCg0rUaNt0aNt9u/+0xmhSyDiQy9iQy9EZMDZXNajZrm1fxoXs2P8aEhnI1NY+fZOHZHxXExPoP4dD1r/u8aa/7vGp6uGtrWtMw116Zmedx1MqJVCCGE/UgCJ+6LRq2y9qUzmXWkZxtIyTSiN5ntHdpdUalUPFDRmwcqevNcx5pcScxgV1Q8u87GcfJaCunZJracjmXL6VhcNCpaVPejY0gA7Wv74+shy3oJIYQoWZLAiSKjUaso567D282FNL2RlAwjepPJZoSrowj282BIKw+GtKpKfFo2e85ZBkEcuZSEwaSw73wC+84noFZBwyAfOoQE0DHEn8o+sqyXEEKI4ldgArdv3z42btzItWvX0Gq1VK9enUGDBlG9evWSik84IJVKhberC96uLuhvTiKcefPLkZpYc/h7uRLWpAphTaqQlm1k//kEdkXFceBCApkGE8euJHPsSjJf/nGO2hU86RASQKeQAGpV8JRlvUQuUq4KIYpCvgncnDlziI+Pp23btsTFxREUFES1atV46aWXGD16ND179izJOIWD0mnV6LSWQRAms0JGtpHUbCPZRsesmfNy1dK1fiBd61uW9Tp8KZFdUXHsiYonKdPAuRvpnLuRzvd7/6Gyjxsdbg6CaFBFlvUSUq4KIYpOvgnc9u3bWbNmDQC9e/fmqaee4j//+Q89evRg2LBhUtCIu6ZRq/B2d8Hb3QWDyYzBaCbLaCJTb0JvMt93QmcwGDgfdZbqNWri5l78TZk6rZq2tfxpW8sf0yMKJ6+msOvmShDXkrO4lpzFikPRrDgUja+7C+1rWwZBtKjuh04rI1rLIilXRUHOnz+Pm5sbVapUsXcowgHkm8CpVCqSkpLw9fUlNjYWs9nSKd3HxwfFEatORKniolHjolHj4aoFT9AbzWTojSRnGu6pmfXYkUOMfmYIsTHX8fH145Mvv+ahrt05dGAfa1evpFLlKgx9ZgTe5XyK4dlYktNGwT40CvZhTJdaXIhLv5nMxRMVm0ZSpoF1x6+z7vh13Fws67l2DAmgbU1/vNykK2pZIeWqyEtWVhaPP/4469atQ61WM3r0aL744gtiY2OJiIggKSmJZ555hsaNG9s7VFGK5PuXY8yYMfTv35+aNWty/vx53n77bQASEhKoV69eScUnyghLU6sOL1ctSRkGUrMNd1UjN+ONKcTGXAcgOSmRN1+dyKyP5zFy6ADrH8l1v61i1YY/ir1fmkqlolYFL2pV8OKZdjW4npzFrijL9CR/RSeTZTCz4+84dvxtWdarWVVfOoQE0CHEnwAvWdbLmUm5KvKyePFi1q1bB4DZbObLL7+kf//+jBkzhvPnzwPw+eefs3fvXpo1a2bPUEUpkm8C16tXL9q3b8/ly5epXr065cqVA6B8+fJ8/PHHJRagKFu0GjUB3q54uWlJyTSQrjcWKpH758I5m+1r0Vf46duvrckbwF9HD3PsyCGaNG9Z1GEXqJKPGwNbBDOwRTBJGXr2nk9g19k4Dv6TgMGkcPCfRA7+k8hnW85Sv7K3dY3WauU9SjROUfykXBV5OXv2bK59q1evtiZvANnZ2SxatIh58+aVZGiiFCuwI46vry+NGjWyFjI5zp07l88VQhQNNxcNgeXcqOzjTjl3lzsOAOjWK8xmu3NoN7y8vXOd5+HpVaRx3i1fDx09G1bivf4NiRzbgbfDHuSR+oF4uVo+S526lsrCnRcYvvhPhi/+k693nuf09RRpXnMiUq6K2/Xr18+mZcDd3Z0OHTrkOs/Ly77llyhd7qnzzXPPPccff/xRxKEIkZubiwY3Fw3l3FxIytCTnm0kr1Rm+swP8fX1Y++uHTRq0ozJr71JzPXrbNmwnpTkJAAeGzCYOnVLTzOVu05D5wcq0PmBChhNZv7vimVZr13n4ohP03MpIYMlBzJYcuAyAV66m3PNBdAk2AetLOvldKRcLbs6d+7MsmXLmD9/Ph4eHrz22mu0b9+er7/+mq1btwJQpUoVxo4da+dIRWmSbwI3c+bMPPcrikJKSkqxBSQsr7HMH2ZLp1UTWM6NLL2JlKzcTatu7u5MeXOGzTW+fuXZuv8of2zeSKUqQbTt0KmEoy48rUZNi+p+tKjux/iuIZy5nnqz31w8lxIyiEvTs/roVVYfvYq3m5a2tfzpEOJPqxrlcXeRZb0cRVkoV6X8ujcDBw5k4MCBNvs2btzIxo0bSUpKok+fPnjn0aogyq58E7hff/2VqVOnotPlXibo999/L9agyqpVq1YxceJErl27xtChQ/nqq69wc3O784VliJtOg5tOg8FkJjXLQEqmgYIGrfqV96f/4KElF2ARUKtU1K9cjvqVy/F8p1pcis+wTk9y+noqqVlGNp2MYdPJGHRaNS1vLuvVrpY/Ph4u9g5fFMCZy9WrV6/yzDPPsHXrVho3bszixYulw/190mg0MrWMyFe+CVyjRo2oU6cOzZs3z3VMOlEWvfj4eJ588kmysrIA+P777wkJCeHNN9+0c2Slk4tGTXlPV7xcXUhI15OhN9o7pGJTzd+DJ/2r8WSbatxIzWbPuTh2nY3j6JVk9EYze87Fs+dcPGoVNA72uTmiNYBK5ST5L22cuVwdP348W7ZsAeD//u//GDp0KKdOnZLaOFGqZOqNpGebKOfu4vDzceabwP33v//F1TXvKQ1y2uRF0Tl69Kg1ecuxf/9+O0XjOHRaNRXLuZKSqSY1y4jBZEYBkhIT0Lq44OXlXE0OFbxd6ds0iL5Ng0jLMrLvQrx1Wa8sg5mjl5M5ejmZ+dvOERLoRaeQADrWCaCGv4f8IS0FnLlc3bdvn832mTNnSEpKws/Pz04ROa7s7GwSExOpVKmSvUNxCkaTZZ7R1CwjeqMZVODt7vjzb+abfmZkZOBeArPZC4vmzZvj4WE7bUSnTqW3z1ZpolKp8PHQEVzegwqeLkydMJpW9WvSql4N5sx6x97hFRsvNy2P1K/I22ENWPVCe2b2a8CjDSrh425pRo2KTWPxnos8991Bnl50gC//OMfx6GSHXI/WWThzuXp7edWgQQNJ3u7Bzz//TOXKlalcuTIdO3YkNjbW3iE5rCy9iRup2UQnZRKXpifbaPmAn+dIOAeUbwI3btw46/fjx48vkWDKMj8/P1asWEH9+vUpV64cL7zwApMmTbJ3WA5n+dIlLP9lCWazGb1ez5efzuHQgX13vtDBubpoaF87gCmP1mXFmHbMHdyEgS2CrM2oV5OyWH7oChN+Ocrgr/by8ca/2Xc+3vJpVJQYZy5X582bx2OPPYaHhwft27dn6dKl9g7J4SQnJzNq1CgSExMB2L17N2+99Zado3IsJrNCaqaB6MRMriVnkpp1b6v7OIJ86xBvnXfq8uXLJRJMWdezZ0/psHqf/vrrr1z7zv99ihat29ohGvvQqFU0qepLk6q+vNClNudupFunJzl/I53EDANr/7rG2r+u4aHT0KZmeTqEBNC6ZnnrfHSieDhzuVqhQgVWr15t7zAc2vnz58nIyLDZl1eZJnLLMphIyzKSrjc6bcJ2uwLXQs3reyFKs549ezJ37lzrtlar5fGwnvh76ki6x3VWHZlKpSIk0IuQQC+Gd6jB1aRMdt8c0Xo8OoUMvYltZ26w7cwNtGoVzav50rFOAO1rB1DeM/dISXF/pFwVBWnYsCHBwcFcuXLFuk8+1OfPbFZI1xtJzTSSbTQ5S8tooeWbwJ0+fZrmzZujKArZ2dnWUVM5c/wcPny4xIIsSy5evMjVq1dp06YNGo3M73W3unXrxsKFC/n8889xc3OjatWqhIaGEhwczKz3P6BB8zak3eU6q86kiq87g1pWZVDLqiRm6Nl7zjII4tA/iRhMCgcuJnLgYiJzN53lwSrl6BASQKeQAIL8nLPfVklz9nI1PT2dQ4cOUb9+fSpUqGDvcByOi4sL69atY8qUKZw/f56mTZuyYsUKvvnmG8LDw5k6daq9QywVzGbFMo1UlgGDqYwW5oBKKUNr9AwYMICVK1faO4x8vfbaa8yePRtFUQgJCWHr1q1UrVrV3mE5rE8++YSXX37Zuu3r68vly5dxcXMnLctIalbBc8iVJRl6IwcuJLIrKo795+NJ15tsjtcM8KRDiD8dQwKoE+iFu05LFV/HSOpK+/u+sEr789izZw99+vQhMTERV1dXvv76a5566il7h+WwoqOjqVWrFnq93rrv559/ZsiQIXaMyr6yDCYy9EbSsowY76PwVgFV/Nxx1Zb+SpKC3vf5DmLYu3ev9fvb+2ps3LixiEITOc6dO2dN3gCioqKYPXu2naNybJs3b7bZTkpK4tChQ7hqNfh7uVLF1wNvNxfusMxqmeCh0/JQ3QpM612flWPbM/vxRoQ1rmxtRr0Ql86P+y4x5sfDDF24n/nbosgymO5wV3E7Zy5Xp0yZYu18n52dzcSJEzEYDHaOynHt2LHDJnmD3GVaWZBlMJGQns2VhAyuJWWSlGG4r+TNmeSbwH344YfW7ydMmGBz7Msvvyy+iMqo6OjoXAuW39oPQty922eB1+l0PPjgg/9ua9VU8P43kZM8zsJFo6ZVjfJM6vYAy0a3Zd7QpgxpVZXgm82osanZ/LT/EhtOXLdzpI7HmcvV2xPShIQEMjMz7RSN42vatGmufWVlZYucJtJrSZaRpEkZBvQ35/gU/8o3gbs1mbg9sShDra4lpl27dtSoUcNm35NPPmmfYJzE1KlTCQsLQ6VSERAQwKJFi/Lsl5OTyFX2cYwq9ZKkVqloUMWH8M61+G5EKxYNb8nIDjUY1CKY0HqB9g7P4ThzuXp7edW7d2/KlStnp2gcX/369fn0008pV64cGo2GZ555hvDwcHuHVayyDCbi07K5nJjBjdRsMg2mMttfuTDuaRSqjJ4qei4uLmzbto0PPviA6Ohohg0bxuDBg+0dlkPz9vbmt99+IzU1FXd3d7TagqfIcNNpqOLiRlKGnuQ7rLFaFqlUKmr4e1LD3xM3Fw3ebrLu6t1y5nJ15syZBAQEsHnzZpo2bcprr71m75Ac3ksvvcQLL7yAwWDA09PT3uEUm4xsIylZRjL1Rqlluwv5/kW7fPkyY8aMyfU9SNNecalRowYLFiywdxhO5caNG6xdu5bg4GC6du16xz+SKpUKP09XPFy1JKYbpEARRcqZy1WNRsPLL79sM3BI3B+TycTGjRtJSkoiLCwMHx8fe4dUpPRGM0kZetKzpZy9F/kmcF988YX1+5EjR9ocu31biNLor7/+omPHjqSkpADwxBNP8MsvvxTqWlethko+GtKzjSRlGMg2Sod9cf+kXBWFZTabeeSRR/jjjz8AqFSpEvv376datWr2DawIKIpCSqahTM7NWZTyTeBat25dknEIUeTmzJljTd4Ali5dyvTp020GMtyJp6sWD52G1CyDjH4S903KVVFYW7dutSZvANevX+fLL7/k/ffft19QRSDLYCIhTU+WfCi+b/kOYti8eTM//fSTdXvQoEF07dqVrl278r///a9EghPifqSnpxdq352oVCrKueuo4utOORmtKu6DlKuisIqq/CoNTGaF1GwDMclZXEvOlOStiOSbwH399deEhoZat/V6PStWrOCHH37g559/LpHghLgfY8aMQa3+91e8devWtGzZ8p7vp9WoCfB2pWI5N1w0ksaJuyflqiisHj16EBISYt12c3Pjueees2NEd8dsVkjLNhCbksWVxAxupGSTrjfKqNIilG8TqsFgoHLlytbtFi1a4Ofnh5+f333P7ZOUlMSkSZOIjo4mKCiITz/9NM/Oma+99hp//PEH/v7+/P7773d9vSjbHnnkEfbs2cOyZcsIDg5m1KhRRTLSz8NVi06rJi5NT4beWASRirKiOMtV4Vzc3NzYu3cvCxcuJCkpiaeffpqGDRvaO6w7ytKbSM02klGGFpW3l3xr4G7tOwQwffp06/cJCQn39aARERG0a9eOjRs30q5dOyIiIvI8b8CAAXz99df3fL0Qbdq04eOPP2bSpEl4e3sX2X21GjWVfNyo4OWKVpZyEIVUnOWqcD4BAQHWJRZLe/KWZTBx/WYTaWqWDE4oCfkmcI0bN2bZsmW59v/yyy80btz4vh50y5Yt9OvXD4B+/frluzxIq1at8qxZK+z1jiYjI4OIiAjeeustjh8/bu9wRCF4u7tY+sa5u+Dg03iJElCc5WppsG3bNl5//XWWL1+O2Wy2dziiBJjNCgnp2VxLziRDpl0qUfk2ob7++uuMGzeONWvW0KBBAwBOnDiBXq9n/vz59/Wg8fHxBAZaZnEPDAy860+e93t9adWjRw927doFwPvvv8/WrVvp2LGjnaMSd6LVqAnwcsVLpyU+XS9Tjoh8FWe5am9fffWVzbx2L774IvPmzbNjRKI46Y1mMvVGUrIMGEySttlDvgmcv78/v/zyC3v37iUqKgqALl260K5du0LdePjw4cTFxeXaP3HixHuL1Mn9+eef1uQNLH1l5s+fLwmcA3HTaaisdSMhQ09qpkE+iYpc7rdcLc3mzp1rs/3VV18xe/ZsPDw87BSRKA7ZRhMpmUbSsg0yIMHO8k3gkpKSAMt6bPXr18+139fXt8Abf/vtt/ke8/f3JzY2lsDAQGJjYylfvnyhAy6K60sjF5fcyxLltU+Ubmq1igAvVzx1GhLTDTJcXti433K1NLu9vNJqtTajwIVj0xvNpGQaSJXErdTIN4EbMGAAKpUKRVG4ceOGtclSURRUKhVbtmy55wcNDQ0lMjKS8PBwIiMj6dq1a4leXxo1bdqU3r17s3btWgA8PT2ZNGmSnaMS98pdp8VdpyUt20ByhlGaVQVQvOWqvb3++usMGzYM5eZf95dffhk3Nzc7RyXul95oJjXLQGqWrA9d2uSbwG3dutX6fb9+/YiMjCyyBw0PD2fixImsWLGCypUr89lnnwEQExPDtGnTWLhwIQCTJ0/mwIEDJCYm0rlzZ8aPH8+gQYPyvd7RRUZGsmbNGq5evUrfvn0JDg62d0jiPnm5uuCp05KWbSQpQy99Rcq44ixX7W3o0KE0bNiQLVu20KxZM7p06WLvkMQ9UhSFDL2JjGwTaXqpcSut8k3gblUUc2fdys/Pj++++y7X/ooVK1qTN4BPPvnkrq53dFqtlv79+9s7DFHEVCoV3m6WRC4xQ09KlhSIoujL1dKgUaNGNGrUiNQsA1kGE65atVM+T2djMivojSayjWYyDSb0RjNmsyL9eEu5QiVwQoj7p1ar8PdyxVNGqwonl5RhwGgyo1arcNGo0WnVuGrU1u/VMnei3WUZTGTojWQZzOiNJhQFSdgcTL4J3OLFi63fx8fH22wDjBgxoviiEsKJuek0VHFxIylDT1Km1MaVJWWnXLXU3pjMCiaziSyD5cOKClCpuJnIaXDVWhI8F40arUYGPBQnvdGM3mTCYDSTZTSTpTdJwubg8k3gbl00d/DgwQ67iK4QpZFKpcLP0xV3nZaENL2MVi0jynq5qgCKAtlGM9lGM6lYkjpUoFGp0N5M5tx0alw1GnRaSeruldFkJstoIttgJlNvwmA2y4dFJ5NvAvfiiy/me1FGRkaxBCNEWePmoqGyrxspmQaSMmX5GWcn5Wpuys1/jIqC8WZtXWqWpaZOq1bhqtXg5qKWmrpCMJkVUrMMpGcb0ZskYXN2Bb4TYmJi+Ouvv9Dr9YClyv+TTz6he/fuJRKcEGWBSqXCx0NHZR93vFy1SO8g5yblauEoChhMCmnZRuLS9FxPzuJKYgbRiRnEpmSRkJ5tGSyhN2E0lb1lu8xmxZLsZhtITM8mJsXy+iSk68k2SvJWFuRbA/ftt9+yYMECqlevjl6v55lnnuGDDz6gb9++rFy5siRjFKLQ/u///o+EhAQ6deqEVutYY3R0WjWB5dwsBXKaHqPUxjkdKVfv3e3Nrzly+tWpVSpctGq0ajUataX2Tq1WoVWr0WkcY+DEjRs3OHjwIC1atLDOEZhDbzSTbbI0iWYbTBhu1rBJKVF25fsXbtmyZfzvf//D19eXq1ev0r17d3788UeaNm1aguEJUXhPP/00P/74IwB169Zlx44duQpBR+Dt6oKbVkNiup70bFkc2plIuVr0chI7s6Jg1JsA2/6ktw6ccNVqcNWWztGwv/32G0888QRZWVnodDqWLFlC3379Sc40kKE3YjQrUqsmbOTbhOrq6mpd1qVKlSrUqFFDCpkSsG/fPn799VdSU1PtHYpD2bdvnzV5Azhz5oxDL6TtorHUxgWUc0Vbiv7IiPvj7OXq1atX+eWXX4j6+4y9Q7FSAPPNmruULAM30rK5lpzJpYR0ohMzuJGaTXKGnky9EaPJbF1JolhjUizNnzlfadkGJk2eTFZWFgB6vZ6JkyYTnZRJcqZlsXhJ3sTt8q2Bu379OjNnzrRux8fH22xPmzateCMrg0aMGGFdQ7ZChQrs3LmTunXr2jcoB3H9+vVc+65evcqrr77KkiVLCA4OZs6cOXTq1MkO0d07b1cX3LWWdVXTsg1SG+fgnLlcXb9+Pf369bP27Zs2czYjwsfaOaq85dUcmzMaVoWlOVatUqHTWmrqNCoVKjWoUKFRq9CqLf/faZLinAlyDSYzBpOCyaxgVhT0RrPNgCUFiLkeY3NtbGwMG9evZe4HM0lJSeaJp4YzbtKrRfo6CMeWbwI3ZcoUm+0GDRoUezBl2cmTJ63JG1j6QsyZM8dmZQqRv27dulG5cmWuXbsGgFqtxtXVlTlz5gCWZC4sLIwrV67g5eVlz1DvmlajpkI5V9yz1SSmy3JcjsyZy9U333zTmrwBfDr7PYYNH4VOp7NjVIWXMxrWUmOnAAp6kxmybc+7PdHTqNWo1bccu3kvk0nBaC58P7X+g4fy4+J/y/sevR9j3MinMBgMAHzy/gyqVq/BYwMG3fuTFE4l3wQuvyWdsrOzbdbzE0UjKSkp177ExMSSD8RBeXp6smvXLj755BMSEhIYOXJkrqXYkpOTOXjwIA899JB9grxPXta+cVIb56icuVy9vbzKzEjHoNc7TAJXWLcnekZz0czhOG3mbKrVqMnB/Xtp1rI1AYEVWbNyuc05e3ZskwROWBVqQh2TycT27duZMmUKDz/8MOvXry/uuJzeoUOHaNq0KRqNhkcffZSaNWvSqFEj63GVSsWoUaPsGKHjqVWrFp9//jlLlizhkUceoXnz5jbHXV1dHb7GI6c2LsDLFY30jXNojlyuZmZm8swzz6DT6ahVqxa///47zz//vM05ffoNxNPBarvtycXFhedeGM+X3y4h/MWJNGnWPNc5DRo1LfnARKlV4DwLf/75J2vWrGH79u00btyYw4cPs2XLFtzd3UsqPqekKApPPPEE586dA2DDhg1MnjyZbdu2MW/ePK5cucJTTz3lsDVFpcV//vMfTpw4werVq6lQoQKffvopFSpUsHdYRcLb3QWdi5r4VFnFwdE4Q7n6wQcf8MMPPwBw4cIFhgwZQnR0NNWrV+d///sftes3YvDTI+0cpWOrXacu02d9xKez3yMjPY2+A59gyDPOstSaKAr5JnCdO3emSpUqDBkyhClTpuDl5UVoaKhDFTKlVWxsrDV5y7F37142b97MV199RWxsLKmpqbRu3RoPDw87Ren4vL29WbVqFenp6bi5uaHRaOwdUpFy1VpWcZA1VR2Hs5Sre/futdlOT09n06ZNzJ8/n927dxNS5wEaNWtJk+Yt7RShc3h21BiefPY5TEYjbg72OyKKX75NqN27dycmJob169ezbds2MjIy7jjiRhROYGAgDzzwgM2+li1b8uyzz3L9+nXMZjPLli3jww8/tFOEzsXT09PpkrccOWuqVvZxx03rnM/RmThLuXr7aG5vb28WL17M7t27AYg6+zeTXhhVIlNyODsXFxdJ3kSe8k3gpk2bxtatWxk+fDj79++nR48eJCQksG7dujK3AHNRU6lULFu2jDZt2uDm5ka/fv0YNmwY2dm2w50OHjxopwiFo8lZU9XXwwUHzAfKDGcpV1999VWef/55vLy8qF+/PitWrODo0aM25/xz4RwpyUl2iU+IsqDAPnAqlYp27drRrl07DAYDO3fuZO3atbzzzjvs37+/pGJ0Sk2aNGHfvn3W7eTkZLy8vEhLS7Pue/jhh+0RmnBQKpWK8p6ueOi0JKbryTRI37jSyBnKVTc3NyIiIoiIiLDue/jhh/npp5+s2/UbNMLH188e4QlRJhR6sUgXFxdCQ0MJDQ21zhYtio6Pjw+RkZG8+uqrREdHM2zYMF566SV7hyUckKU2zp3UbANJMm9cqeZM5ep///tfDAYDmzdvpl6DRsz4cK69QxLCqd0xgZs3bx7jx4+3bs+dOxdPT08GDRqEn598uipKXbt25fDhw/YOQzgJb1cXPF20JGfqScky2sz8LuzLGcvV8uXLs3TpUgAuJ6TLBwchitkd54G7fd6sRo0aodVqef/994stKCFE0VCr/x3k4KErdIW7KGZSrgoh7tcdS/TQ0FCb7UceeaTYginrsrKy+Pnnn4mOjmbgwIHUq1fP3iEJJ6HTqqnk40ZqpoHEDD1GqY2zK2ctV3ft2sXmzZsJrl2Ph3v0dsgRtkI4inxr4BISEmy2V69ezcyZM1m6dKkMDS8mjz76KCNHjuTNN9+kadOmueZaEndHURQ+/vhjWrRoQd++ffnrr7/sHZLdebu7UNnXHS9XLfKnteQ5c7n6zTff0KlTJ9555x2ef2Yo701/zd4hObyD+/cy/In+PN7zYZb99J29wxGlTL4J3HPPPWf9/osvvuC3336jQYMG7N69W6r5i8HBgwfZvn27dTs7O5t58+bZMSLHFxERwSuvvMLhw4f57bff6Natm8N3FC8KLho1geXcqFDOFZ2mUKvpiSLizOXqxx9/bLP90+KFZGZk2Ckaxxcfd4MRT/Rn57bNHD10kNcmvciWDevsHZYoRfJtQr310+CmTZv46aef8PDwoE+fPgwYMKBEgitL1Orcf0iddfLZkvLbb7/ZbMfExHDgwAE6d+5sp4hKFy9XFzx1WlKzDCRlGKRZtQQ4c7l6exmmUqulCfU+7Nm5nYwM27kBN29YR9cevewU0f0xKwoGoxmDSUFvMmO4+aW/uc9gMlv3642K9XjOtt7m/H+vydmnv3XfbY9z6zV6kxmAiY/UYexDIXZ+Ve5PvglcVlYWJ0+exGw2YzKZrEs6ubi45JlsiPvTvHlzevTowYYNGwBwd3eXaUTuU7169Vi37t9PrFqtlpAQx37DFjWVSkU5dx1eri4kZepJliW5ipUzl6tTp07lmWeesSapI8LHygoC96F2yAOF2pcfRVEwmnMSmluSmVsSnPySptsTnkIlTQU8jt5kLnWj4E9dS7F3CPct3wSuQoUK1ip9Hx8fYmNjCQwMJDExUWqGismaNWtYuXIlV69epX///tSoUcPeITm01157jT179rBv3z7c3d354IMPqFKlir3DKpXU6n8nAU5I05NllEmAi4Mzl6tPPfUUDz74IJs3byaodn3adS57E5Hn1DLpb0l+DMZbE5pbkybFJlm6PeExmDxp1e85Dq35DrPJSFDDdpwLaMfUX4+hvzVpyus+JkscpStlujMV4KJV46JRodOocdGo0Wkt/7toVDf/V6O7+f2/x/49rtOqb16runmvW665eczdRUO3BhXt/XTvW74J3A8//JDn/nLlytnMti2KjouLC0888YS9w3AaAQEBLFiwgIULF1K7dm1GjRpl75BKvZwluZIy9CRJbVyRc/ZytXnz5jRv3rzE5oG7tZZJf3vz3K2JlLUGyWxJfmy286o5Um47bi5U0lTUtUzm6t3wbq/HnJmM0qw3+y9nAEXbr1CjVv2bMFmTn5v/a1U2CVJ+SZXOui/vxEp3y31vTax0tz2OVq0qkWZ3FZa+wI7urieG0mg0uEu1uHAA27dv55FHHsFoNAKwcuVKdu7caeeoSj+VyjJ3nLvUxpUYZypXz1xP5eTVJDL0tyVNNk1xtonVrbVGt15z+/l5JU2OpjC1TDqtGo1iYtN7L5Jy7SIAmX9tZOjMH6lcIyRXzZJtwpRP0nRLgnZrDZVa+ik6rHua2bN///6sWrWqqGMRokjNnz/fmryBZY6qQ4cO0aJFCztG5Tjyqo1LSU5i0YL5XL54nqGDH2fw4MH2DtNpOEO5uvLwFSYv+z97h5GLVq2yqf35N4n5t/bn1u2Cj+fVhHd7TdRtTX63NOtpClnLtGXDOn69mbwBGLOz0JzdyuhnexTjK+X8fo/8lU3rf6dJw/q8+vLL+Pj42Duke3ZPCZyjFzKibNDpdIXaJ/J3e23c4GGDOHRgHwCRvy4jJSVFmqaLiDOUq5V83PB205KebSygb9K/NUQF1RjdXquk097aB8o2sSooadJqVA5Zy+TikruscnFxsUMkzmPJd4t481XL4MDfV8HO7dttpu9yNIVO4NLS0rh48SJVq1Z16IxVlB2TJ08mMjKS9HTLUPy+ffvSqFEjO0flmNxcNKTHXbEmbzkWL14sCdx9cLZytX3tAI691Z0riRkO2bxZmnTo8jDNWrbmyMEDAJT392fYcHmv3Y+Vv/xos71jxw7Onz9PrVq17BTR/cm3F98rr7xinTV8586d9O7dmzlz5tCvXz/Wr19fYgEKca+aN2/OmTNnmD9/PqtXr+bXX3+1d0gOzc/PD63W9jNfYGCgnaJxTGWhXJW534qGRqNhSeR65n39PTPnfMaGXQepWr2GvcNyaOUDAmy2dTodvr6+9gmmCORbA3fmzBnKly8PWPoS/fTTTwQHB5OQkMDw4cPp2bNniQUpxL0KCgpi7Nix9g7DKQQEBPDGG2/wzjvvAJaEbvr06XaOyrFIuSruhk6no9dj/e0dhtMY/8pr/LlvLynJSQBMmzbN+n50RPkmcGazmbS0NLy8vFCpVNb5s8qXL4/JJKPShCiL3n77bYYMGUJUVBRdunTB29vb3iE5FClXhbCfRk2asePwCf7cu5s2TRvQqMGD9g7pvuSbwI0bN45nnnmGJ598kubNm/PSSy/RtWtX9u3bR6dOnUoyRiFEKVKvXj3q1atn7zAckpSrQtiXt3c5unbvSRU/x5+2J98ErlevXjRo0IBly5Zx8eJFTCYTR44coXfv3lLQCCHEPZByVQhRVAochVq9enVeffXVkopFiCK3atUqlixZQnBwMK+++qospSXsTspVUVgXz59j0YLPSUlJZtCTz9Ch80P2DkmUIvkmcO+//z7du3eXSU+Fw1qxYgWDBg2ybq9bt46TJ086/JqTwnFJuSoKKy0tlcF9uhEfdwOAtZG/8nPkelq2bW/nyERpkW8Ct3r1av78808SExPp2bMnffr04cEHHbvDnyhbbl938u+//2b//v20by8FoLAPKVdFYe3ctsWavIFlAMzqX5dJAies8k3gKlWqxMqVK7l48SJr167l1VdfxWQy0adPH3r37k3NmjVLMk4h7lqlSpVstlUqFRUrVrRTNEJIuSoKr0KF3HMsVgiU8kv8K9+JfHMmY6xRowbjxo1j7dq1fPrpp2RnZxMeHn5fD5qUlMSIESPo3r07I0aMIDk5Oc/zXnvtNdq1a0efPn1s9s+bN49OnTrRt29f+vbt69BLYYjiM3XqVKpWrWrdnjRpErVr17ZjRKKsK85yVTiXlm3bEzbg3y4gderW46mRz9sxIlHa5FsDpyi5l0HJmT7g5Zdfvq8HjYiIoF27doSHhxMREUFERESenXoHDBjAU089xX/+859cx4YPH85zzz13X3EI51azZk2ioqLYuXMnwcHB1K1b194hiTKuOMtV4Xw+XbCI0S9OIiUlmZZt2kn/XWEj3xq4n376qdgedMuWLfTr1w+Afv36sXnz5jzPa9WqlVOsDyjsR6fT0bVrV0neRKlQnOWqcE71GzaiTfuOkryJXPKtgfP09ERRFI4dO0ZMTAwqlYrAwEAaN25832vdxcfHW9dQDAwMtK4NeDd++uknIiMjadiwIVOnTpVETwhR6hVnuSqEKFvyTeB27drFO++8Q/Xq1a0dv69fv86lS5d466236NixY4E3Hj58OHFxcbn2T5w48f4iBoYOHcrYsWNRqVR89tlnfPDBB7z//vv3fV8hhChO91uuCiFEjnwTuPfee4/FixcTHBxss//y5cuEh4ezfv36Am/87bff5nvM39+f2NhYAgMDiY2NvevFZAMCAqzfDxo0iDFjxtzV9UIIYQ/3W64Kx6ECVCpQq1Ro1CrMioKiWPaZlX/7Qyo3z82hKJZ9QtxJvgmcyWTKNQ0DQMWKFTEajff1oKGhoURGRhIeHk5kZCRdu3a9q+tzkj+AzZs3U6dOnfuKRwghSkJxlquieN2akKlVKtRqS2Kmtu4DkwIooNOq0Gk0uGjVaNR5N42bzQoKlkROfbP53KwoGM0KJrOC3mjCYFIwmBTMitma9EmCJ3Lkm8A9/vjjDBw4kF69elG5cmUArl27xrp16xg4cOB9PWh4eDgTJ05kxYoVVK5cmc8++wyAmJgYpk2bxsKFCwGYPHkyBw4cIDExkc6dOzN+/HgGDRrERx99xOnTpwEICgpixowZ9xWPEEKUhOIsV8X9UQEatQo3Fw0aterf2rObyZpWo0KnURdZX0W1NbH7935qVGhvjlXwdM3959loMmM0K5jNCkazGb1RwWA2YzCaMd1M7kTZoVLyGtd+07lz59iyZQsxMTEoikKlSpUIDQ0lJCSkJGMsMgMGDGDlypX2DkMIUYJK2/v+XsvV0vY8CnI5IR2DyTGyCZUKPHVavFy1eOSRNDmKbKMJo8mS3JkVBbMCppuJntFk+V8BkBo8VEAVP3dctaV/ZG9B7/sCf1tr164tE58KIUQRknLV/lQq0Gk0eOjUeLq6oNPmO6OWw3DVaigo/1RymmdNCgbFjN5gJstgRm8ySVLnoPL9rd2xY4f1+9TUVN544w3CwsJ4+eWX8xxdKoQQomBSrtqPRq3C01VLBW9Xgv08CPJzx8/T1SmSt8JQqVS4aNS46TR4u7rg7+VKkJ87VXzdqeTjRsVybgR46Sjn5oJOoyafrnuiFMn3N3fu3LnW7z/44AMCAgJYsGABjRo1Yvr06SUSnBBCOBMpV0uOCtCqVbhq1fh5uBDk607Fcm54u7ngoikbSVthuGo1uOu0eLpqKeeuI8DbleDyHgT7eVCpnBvlPXV46LRob/YLFKVHoRr8jx8/zurVqwHL/G6rVq0q1qCEEMLZSblaPFSAu05LOTfH7tNmb1qNGq1Gjcct+7KNJvQmS/Or3mS2DqqQwRP2ke/HkPj4eBYvXsyiRYtIS0uzWcPPbDaXSHBCFFZsbCyPPfYYHh4edOjQgRMnTliPnTx5UpqnRKkg5WrxUavAy1VLZR9Lk6CjJW+zZs0iMDCQ4OBgIiIirPuvX7/O2bNn7RjZv1y1/za/VvZxp2p5T6rerKnz83DBXadBp1FbR/GK4pVvAjd48GDS09PJyMigf//+JCYmAnDjxg3q169fYgEKURjjx49nzZo1ZGZmsmfPHoYMGUJMTAzNmjWjQYMGVKlShQ8//NDeYYoyTsrVoqXCMoK0gpcrQX4eBJZzw01X+kcW3u7333/njTfe4MaNG0RHRzN69GiOHDnCK6+8QlBQEA888AAPPfQQqamp9g41F61GjYerFj9PS1IXXN6DauU9CPK1JHb+npZ+dbdOzyKKRr4fUV588cU891eoUEH+EIpSZ+fOnTbbx48f5+233+bo0aMAGAwGXn/9dZ588slcs+ALUVKkXC0aKiyjSP08XRyupi0vt5dfAEuWLOHjjz+2bm/fvp0vv/ySKVOmlGRo90SlUlkmM85jgIjemNP0apnHTm8yYzCZMd+coVhaYwuvwN/8nPmKYmNjAcvC8127dpUh8KLUadu2rU0fonr16nH58mWbc0wmExcuXJAETtiVlKt3L2cVBI1ahbuLFg+dxikStxxt27bNtc/Pzy/XvqioqJIIp1jptOo7J3YmheybU5xI/7r85duEGhERweTJkwFo1KgRjRo1AiyrI9zaPi9EafD555/zyCOPoFKpaNq0KT///DMDBgywOScoKIg2bdrYKUIhpFwtjJzRo56uWnw9XKhQzpXKvjf7W5X3JMDb1amSN4D+/fvzxhtv4O3tjb+/P59++imjR4/G29s713nOSqe1NMWWc9cRcHOKkyBfDwK8LKNgpfk1t3zfBb/++iu///47Li4uNvuHDx9Onz59CA8PL/bghCisKlWqsGnTJhRFsS5107RpU7KysliyZAnBwcG8/fbb6HQ6O0cqyjIpV/OWU8Pm5eqCl5sWV23RLVnlKGbOnMnMmTNt9m3dupX33nuPxMREwsPD6dmzp52isw9LbZ2Ocu6W7ZxaOoPJjMGkkG201NiV1ebXfBM4lUpFbGwsQUFBNvtv3LhR5t5YwnHc/rs5duxYxo4da6donE9mZiaLFi3i7Nmz9O3bl4cfftjeITkUKVf/ZVkNQY2HToO7ixadVn3L+qACoGXLljK9zC3yan5VFAWDScFkNmNUFAzGmytMGC3Nr7cndXt2bmfL/9bSrNGDhD8/Cjc3t5J7AkUs3wTu9ddfZ/jw4VSvXt266PLVq1e5dOkSb775ZokFKIQoPQYOHMi6desA+Oyzz1i6dCmDBw+2c1SOo6yWqyosi7dr1Sp0Wg06jQpXFw1uLo43YlSULjkDJqw9wlwt/5nNClkGExkGExnZRkxmhdUrlzNpzEgAvgU2bdzAmjVr7BF2kcg3gevcuTMbNmzg2LFjNosuN2rUCI1G3nRClDUXL160Jm855s+fLwncXShr5WrO8lWeOkuyVtZqGYX9qNUqPFwtkzkrnjqyDCZ++e5rm3N+//13Ll26RLVq1ewU5f0psCeoWq0mODgYFxcXVCoVgYGBTlnICCHuzM3NDbVabTPhrKenpx0jckxloVzVqNV46NR4uznHQvHCsalUKtx1Wny8vWz2azQa52xCPXXqFG+99RapqalUqlQJRVG4fv065cqV46233qJBgwYlGWeZkJKSwtdff010dDRDhgyhVatW9g5JCKtKlSoxduxYPv/8cwDc3d15/fXX7RyVY3H2cnXNmjVs3ryZJk2a8Oyzz6KRNUdFKfLGG2+wY8cOsrKyAMu8jIGBgXaO6j4o+XjssceUo0eP5tp/5MgRJSwsLL/LSrX+/fvbO4R8mc1mpWXLlgqWPpeKRqNRtmzZYu+whMhlx44dyjfffKNcuXLF3qEUSml6399PuVqankdePv30U2v5BSijRo2yd0hC5HL58mXlm2++UXbu3GnvUAqloPd9vh+PMjMzadKkSa79TZs2JTMzs7jyyTJr//79HDx40LptMplYsGCBHSMSIm+dOnVi5MiRuUZSijtz5nJ1/vz5Ntvffvst6enpdopGiLwFBwczcuRIOnbsaO9Q7luBgxjCw8Pp168flSpVAiyL6kZGRtKpU6cSC7CsyKsvkfQvEsK5OHO56uHhYbPt6uqKVutcE+4KUZrk++6aNm0a27dvty75oigKFStWZNiwYXTp0qUkYywTGjVqxMCBA1mxYgUAPj4+vPzyy3aOSghRlJy5XH3rrbcYNGgQJpMJsPQ3cnV1tXNUQjivfBO47OxsGjVqlKtQiY+PJzs7W96YxWDZsmVs3ryZ6Oho+vTpQ0BAgL1DEkIUIWcuV/v378+ZM2fYtm0bTZs2pWXLlvYOSQinlm8fuJkzZ9r0ycqxe/duZs2aVaxBlVUqlYpu3boxfPhwSd6EcELOXq7Wrl2bUaNGSfImRAnIN4E7dOgQ3bt3z7X/sccey7MAEkIIUTApV4UQRSXfBE5R8l8W9taJPIUQQhSOlKtCiKKSbwLn7+/PsWPHcu0/duwY5cuXL9aghBDCGUm5KoQoKvkOYpgyZQoTJ06kf//+1tnBjx8/TmRkJHPnzi2xAIUQwllIuSqEKCr51sA1btyYZcuWoSgKq1atIjIyEoDZs2dbvxdCCFF4Uq4KIYpKgbMsBgQEMGHCBE6ePMnvv/9OZGQkf/75Jz169Cip+IQQwqlIuSqEKAr5JnAXLlxg7dq1rF27Fl9fX3r16oWiKPzwww8lGZ8QQjgNKVeFEEUl3wSuZ8+etGzZkgULFlC9enXAsradEEKIeyPlqhCiqOTbB27evHkEBATwzDPPMG3aNPbu3VvgEHghhBAFk3JVCFFU8q2B69atG926dSMjI4PNmzfz7bffEh8fz1tvvUW3bt3o2LFjScYphBAOT8pVIURRybcGLoeHhwePPfYYX331Fdu3b6d+/fpERESURGxCCOGUpFwVQtyvAkeh3s7X15chQ4YwZMiQ4opHCCHKFClXhRD34o41cEIIIYQQonSRBE4IIYQQwsFIAieEEEII4WAkgRNCCCGEcDCSwAkhhBBCOBhJ4IQQQgghHIwkcEIIIYQQDsYuCVxSUhIjRoyge/fujBgxguTk5FznXLt2jaeffpqePXvSu3dvvvvuu7u6XgghhBDCWdklgYuIiKBdu3Zs3LiRdu3a5TkDuUajYerUqaxfv56lS5eyZMkSoqKiCn29EEIIIYSzsksCt2XLFvr16wdAv3792Lx5c65zAgMDadCgAQBeXl7UqlWLmJiYQl8vhBBCCOGs7JLAxcfHExgYCFgStYSEhALPv3LlCqdOnaJJkyb3dL0QQgghhDO5q7VQ78bw4cOJi4vLtX/ixIl3dZ/09HQmTJjA66+/jpeXVxFFJ4QQQgjhuIotgfv222/zPebv709sbCyBgYHExsZSvnz5PM8zGAxMmDCBsLAwunfvftfXCyGEEEI4I7s0oYaGhhIZGQlAZGQkXbt2zXWOoii88cYb1KpVixEjRtz19UIIIYQQzsouCVx4eDi7d++me/fu7N69m/DwcABiYmJ4/vnnATh06BCrV69m37599O3bl759+7J9+/YCrxdCCCGEKAuKrQm1IH5+fjbzuuWoWLEiCxcuBKBly5acOXPmrq4XQgghhCgLZCUGIYQQQggHIwmcEEIIIYSDkQROCCGEEMLBSAInhBBCCOFgJIETQgghhHAwksAJIYQQQjgYSeCEEEIIIRyMJHBCCCGEEA5GErhSZPv27bRu3ZqgoCCmTJmC0Wi0d0hCCFEoiYmJDBs2jIoVK9KjRw+ioqLsHZIQTs0uKzGI3FJSUnjsscdISUkB4KOPPqJy5cpMmjTJzpEJIcSdvfTSSyxZsgSAjRs3MnjwYA4fPmznqIRwXlIDV0ocPHjQmrzl2LJli52iEUKIu3N7eXXkyBESExPtFI0Qzk8SuFKiYcOG6HQ6m33Nmze3UzRCCHF3WrRoYbNdq1YtfH197ROMEGWAJHClRGBgIIsWLaJChQqo1WoGDBjAf/7zH3uHJYQQhfLf//6XNm3aABASEsKPP/6ISqWyc1RCOC/pA1eKDBs2jCFDhpCdnY2Hh4e9wxFCiEKrUaMG+/btIzU1FS8vL0nehChmksCVMhqNRpI3IYTD8vb2tncIQpQJ0oQqhBBCCOFgJIETQgghhHAwksAJIYQQQjgYSeCEEEIIIRyMJHBCCCGEEA5GEjghhBBCCAcjCZwQQgghhIORBE4IIYQQwsFIAieEEEII4WAkgRNCCCGEcDCSwAkhhBBCOBhJ4IQQQgghHIwkcEIIIYQQDkYSOCGEEEIIByMJnBBCCCGEg5EETgghhBDCwUgCJ4QQQgjhYCSBE0IIIYRwMJLACSGEEEI4GEnghBBCCCEcjCRwQgghhBAORhI4IYQQQggHIwmcEEIIIYSDkQROCCGEEMLBSAInhBBCCOFgtPZ40KSkJCZNmkR0dDRBQUF8+umn+Pj42Jxz7do1pkyZQlxcHGq1msGDB/Pss88CMG/ePJYtW0b58uUBmDx5Ml26dCnx5yGEEEIIYQ92qYGLiIigXbt2bNy4kXbt2hEREZHrHI1Gw9SpU1m/fj1Lly5lyZIlREVFWY8PHz6c1atXs3r1aqdK3g4cOMCqVatIS0uzdyhCCHFXrl+/zvLly/n777/tHYoQTs8uCdyWLVvo168fAP369WPz5s25zgkMDKRBgwYAeHl5UatWLWJiYkoyzBL33HPP0aZNGwYMGEDt2rWlEBRCOIwNGzZQo0YNBg8eTL169Zg3b569QxLCqdklgYuPjycwMBCwJGoJCQkFnn/lyhVOnTpFkyZNrPt++uknwsLCeO2110hOTi7WeEvCqVOnWLRokXU7NjaWjz76yI4ROY+srCzMZrO9wxDCqb3xxhtkZ2cDoCgK06ZNQ6/X2zkqx2cymeR1FHkqtgRu+PDh9OnTJ9dXXrVtBUlPT2fChAm8/vrreHl5ATB06FA2bdrE6tWrCQwM5IMPPiiOp1Ci8kpiExMT7RCJ80hPT2fQoEF4enoSFBTEsmXL7B2SEE7r9vIqLS1NEo/79OWXX1KhQgW8vLx4/vnnMRqN9g5JlCLFNojh22+/zfeYv78/sbGxBAYGEhsbax2McDuDwcCECRMICwuje/fu1v0BAQHW7wcNGsSYMWOKLG57adeuHQ0aNODEiRMAqFQqRo4caeeoHNvs2bNZsWIFYOmb88wzz9C1a1f8/f3tHJkQzmfUqFG8/vrr1u0hQ4ZYP3SLu/f3338zbtw4FEUB4Ouvv6Z58+a88MILdo5MlBZ2GYUaGhpKZGQk4eHhREZG0rVr11znKIrCG2+8Qa1atRgxYoTNsZzkD2Dz5s3UqVOnROIuTmq1mj/++IN58+YRHR3Nk08+SWhoqL3DcmgHDx602c7Ozuavv/7ioYcesk9AQjix1157japVq7J582aaNm0qicZ9OnTokDV5y3F7mSbKNrskcOHh4UycOJEVK1ZQuXJlPvvsMwBiYmKYNm0aCxcu5NChQ6xevZoHHniAvn37Av9OF/LRRx9x+vRpAIKCgpgxY4Y9nkaRCwgI4J133rF3GE7j4YcfZv369dbtcuXK0aJFCztGJIRze+qpp3jqqafsHYZT6NixI1qt1qbZ9OGHH7ZjRKK0sUsC5+fnx3fffZdrf8WKFVm4cCEALVu25MyZM3leL537RWFMnDiRq1evsmTJEoKDg5kzZw7e3t72DksIIe6oatWqLF++nGnTppGYmEh4eDjDhg2zd1iiFLFLAidESXBxcWHu3LnMnTvX3qEIIcRd69evn3XKLSFuJ0tpCSGEEEI4GEnghBBCCCEcjCRwQgghhBAORhI4IYQQQggHIwmcEEIIIYSDkQROCCGEEMLBSAInhBBCCOFgJIETQgghhHAwksAJIYQQQjiYMrUSQ3R0NAMGDLB3GEKIEhQdHW3vEIqElF9ClD0FlV8qRVGUEoxFCCGEEELcJ2lCFUIIIYRwMJLACSGEEEI4GEnghBBCCCEcjCRwQgghhBAORhI4IYQQQggHU6amESltEhMTGT58OABxcXGo1Wri4uJ44IEHMBgMxMXF4eXlhbe3N35+fnz77bd2jdcR5PWali9fHoA+ffrw66+/otVqUavVjBw5kn79+tkvWAdx48YNZs2axV9//YVOpyMoKIjXX38dgFmzZnHx4kW0Wi0PPPAAb775JgEBAXaOWJQUKcOKlpRfRc+pyy9FlAr//e9/la+//tpm33/+8x9l/fr1dorI8d36mi5ZskQZOXKkkpqaqiiKoqSkpCgrV660Z3gOwWw2K4MHD1aWLFli3Xfy5Enlzz//VLp166Zs2bLFun/v3r3KmTNn7BGmKAWkDCtaUn7dP2cvv6QGTpQJX331Fd9//z1eXl4AeHt7079/fztHVfrt27cPrVbL0KFDrfvq16/PihUraNq0KaGhodb9bdu2tUeIQjg9Kb/ujbOXX9IHTji9tLQ00tPTqVatmr1DcThnz56lQYMGhd4vhChaUn7dO2cvvySBE2WCSqWydwhCCHFPpPwSeZEETjg9Ly8v3N3duXz5sr1DcTh16tThxIkTufaHhITkuV8IUbSk/Lp3zl5+SQInyoTw8HDeeecd0tLSAEuzxNKlS+0cVenXtm1b9Ho9y5Yts+47duwY1atX58iRI/zxxx/W/Tt27ODMmTN2iFII5ybl171x9vJLEjhRJjz55JO0adOGxx9/nD59+vDUU0/h5uZm77BKPZVKxeeff87u3bt55JFH6N27N59//jmBgYEsWLCAH374ge7du9OrVy9WrVqFv7+/vUMWwulI+XVvnL38UimKotg7CCGEEEIIUXhSAyeEEEII4WAkgRNCCCGEcDCSwAkhhBBCOBhJ4IQQQgghHIwkcEIIIYQQDkYSOOFQZs+ezaOPPkpYWBjjxo0jJSUFgN27dzNgwADCwsIYMGAAe/futV4zd+5cunTpQrNmzWzupdfrmThxIt26dWPQoEFcuXLF5nhaWhqdOnVixowZ1n1Tp04lNDSUvn370rdvX06dOlWMz1YI4Uyk/BJFSRI4USopioLZbM61v0OHDvz++++sWbOGGjVq8NVXXwHg5+fHl19+yZo1a/jggw+YMmWK9ZqHH36Y5cuX57rX8uXLKVeuHJs2bWL48OHMmTPH5vinn35K69atc103ZcoUVq9ezerVq6lfv/79PlUhhJOR8kuUBK29AxAix5UrV3j++edp06YNR48exdvbm8TERFQqFY8//jjDhw+nY8eO1vObNm3K//73PwAefPBB6/46deqg1+vR6/XodDqaNm2a5+Nt3bqVF198EYAePXowY8YMFEVBpVJx/Phx4uPj6dSpE8ePHy++Jy2EcApSfomSJjVwolS5cOEC/fr1Y+bMmWi1Wuun1QEDBuQ699dff6Vz58659m/YsIH69euj0+kKfKyYmBgqV64MgFartRa4ZrOZ2bNn23wKvtXcuXMJCwtj1qxZ6PX6e3iWQghnJOWXKEmSwIlSpUqVKjRt2pSqVaty+fJl3n33XXbs2IGXl5fNeV9++SUajYbHHnvMZv/Zs2eZM2eOTb+P/OS1CIlKpWLJkiV07tzZWjjeavLkyfzvf//j119/JTk5mYiIiLt8hkIIZyXllyhJ0oQqShUPDw8AfHx8WL16Nbt27WLJkiWsX7+e999/H4BVq1bxxx9/8O2336JSqazXXr9+nRdffJHZs2dTrVq1Oz5WpUqVuHbtGpUqVcJoNJKamoqvry9Hjhzh0KFD/Pzzz6Snp2MwGPDw8OCVV14hMDAQAJ1Ox4ABA1i0aFExvApCCEck5ZcoSZLAiVIpISEBnU5Hjx49qFatGlOnTgVgx44dLFy4kB9//BF3d3fr+SkpKYSHhzN58mRatGhRqMcIDQ1l1apVNGvWjA0bNtC2bVtUKhUff/yx9ZyVK1dy/PhxXnnlFQBiY2MJDAxEURQ2b95MnTp1ivBZCyGcgZRfoiRIAidKpdjYWF577TXrSK7JkycD8O6776LX6xkxYgQATZo0YcaMGfz4449cunSJL774gi+++AKARYsW4e/vz4cffsjvv/9OZmYmnTt3ZtCgQYwfP56BAwfy6quv0q1bN3x8fJg7d+4d43rllVdITExEURTq1avHO++8U0yvgBDCUUn5JUqCSsmrIV0IIYQQQpRaMohBCCGEEMLBSAInhBBCCOFgJIETQgghhHAwksAJIYQQQjgYSeCEEEIIIRyMJHBCCCGEEA5GEjghhBBCCAcjCZwQQgghhIP5f+4J/7b8bq4xAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\n",
    "ax1, ax2 = axes\n",
    "\n",
    "im_coef, im_p = spearmanr(concat_imputed_df['genotype'], concat_imputed_df[genepair])\n",
    "# sns.violinplot(x=concat_imputed_df['genotype'],  \n",
    "#                y=concat_imputed_df[genepair], \n",
    "#                ax=ax1,\n",
    "#               inner=None)\n",
    "sns.swarmplot(x=concat_imputed_df['genotype'],  \n",
    "               y=concat_imputed_df[genepair], \n",
    "               ax=ax1,\n",
    "              color='black')\n",
    "sns.regplot(x=concat_imputed_df['genotype'],  \n",
    "               y=concat_imputed_df[genepair], \n",
    "               ax=ax1, scatter=False)\n",
    "ax1.set_title(f'Imputed r={im_coef:.2f}; pvalue {im_p:.4f}')\n",
    "ax1.set_xticklabels([f'{refallele}{refallele}', \n",
    "                     f'{refallele}{altallele}',\n",
    "                     f'{altallele}{altallele}'])\n",
    "ax1.set_xlabel(snp_id)\n",
    "\n",
    "coef, p = spearmanr(concat_df['genotype'], concat_df[genepair])\n",
    "# sns.violinplot(x=concat_df['genotype'],  \n",
    "#                y=concat_df[genepair], \n",
    "#                ax=ax2,\n",
    "#               inner=None)\n",
    "sns.swarmplot(x=concat_df['genotype'],  \n",
    "               y=concat_df[genepair], \n",
    "               ax=ax2,\n",
    "              color='black')\n",
    "sns.regplot(x=concat_df['genotype'],  \n",
    "               y=concat_df[genepair], \n",
    "               ax=ax2, scatter=False)\n",
    "ax2.set_xlabel('')\n",
    "ax2.set_title(f'Not Imputed r={coef:.2f}; pvalue {p:.4f}')\n",
    "ax2.set_xticklabels([f'{refallele}{refallele}', \n",
    "                     f'{refallele}{altallele}',\n",
    "                     f'{altallele}{altallele}'])\n",
    "ax2.set_xlabel(snp_id)\n",
    "plt.savefig(example_savedir/f'{snp_name}_ref{refallele}_alt{altallele}_{gene1}_{gene2}.{celltype}_{datasetname}.full.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PosixPath('/groups/umcg-lld/tmp01/projects/1MCellRNAseq/GRN_reconstruction/ongoing/coeqtl_mapping/output/examples/rs221045_C_refT_altC_AC005076.5_ARHGEF19.monocyte_onemillionv2.full.pdf')"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "example_savedir/f'{snp_name}_ref{refallele}_alt{altallele}_{gene1}_{gene2}.{celltype}_{datasetname}.full.pdf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
